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Abstract 

Research  into  improved  calibration  targets  for  measurement  of  radar  cross-section 
has  created  a  need  for  the  ability  to  accurately  compute  the  scattering  from  perfectly 
conducting  bodies  of  revolution.  Common  computational  techniques  use  moment-method 
codes  that  employ  subdomain  basis  functions  to  expand  the  unknown  current  density  on 
segments  of  the  surface.  This  approach  has  its  shortcomings.  Large  numbers  of  basis  func¬ 
tions  are  required  to  accurately  describe  the  behavior  of  the  induced  current  density.  Also, 
increasing  the  number  of  basis  functions  to  improve  accuracy  after  an  initial  computation 
requires  re-computation  of  previous  results  and  lost  processing  time. 

This  research  involves  using  basis  functions  to  expand  the  unknown  currents  that 
have  as  their  domain  the  entire  length  of  the  surface.  Entire-domain  basis  functions  are 
better  able  to  model  the  current  density  on  a  smooth  surface,  and  functions  can  be  chosen 
that  closely  model  the  functionality  of  the  expected  current.  Fewer  modes  are  required 
resulting  in  smaller  matrix  sizes.  In  addition,  accuracy  can  be  increased  by  incrementally 
adding  entire-domain  modes  while  retaining  previously  computed  results  saving  significant 
computation  time. 

Electric-field  integral  equations  are  developed  and  solved  by  an  entire-domain  im¬ 
plementation  of  the  Method  of  Moments  for  a  perfectly  conducting  sphere.  Comparison 
is  made  to  the  exact  Mie  series.  Convergence  in  fewer  modes  is  demonstrated  over  an 
equivalent  application  of  subdomain  pulses.  Matrix  fill  time  saves  as  much  as  hours  over 
subdomain  discretization. 
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COMPUTATION  OF  SCATTERING  FROM  BODIES  OF  REVOLUTION 
USING  AN  ENTIRE-DOMAIN  BASIS  IMPLEMENTATION 
OF  THE  MOMENT  METHOD 


L  Introduction 


Throughout  history,  many  military  operations,  large  and  small,  have  been  successful 
due  to  efforts  to  ensure  surprise  and  improve  survivability.  In  modern  times,  this  has  not 
changed.  Stealth  technology  is  essential  in  both  of  these  areas.  Incorporating  stealth  on 
today’s  complicated,  hi-tech  battlefield,  however,  requires  technically  advanced  and  costly 
equipment  and  weapons.  Although  achieving  military  stealth  involves  many  scientific 
and  non-scientific  disciplines,  one  area  of  intensive  research  centers  on  reducing  the  radar 
cross-section  (RCS)  of  aircraft  and  other  weapon  systems.  In  the  process  of  design,  test, 
modification,  and  use  of  these  systems,  measurement  of  the  systems  or  components  of 
the  systems  on  controlled  ranges  is  required.  The  RCS  must  be  determined  with  extreme 
accuracy,  ensuring  that  the  RCS  measured  is  that  of  the  system  and  not  of  the  system 
interacting  with  the  background  of  the  measurement  range.  Any  effects  of  the  measurement 
chamber  or  background  must  be  subtracted  out.  Calibration  involves  a  transfer  process 
whereby  the  corrections  for  the  radar  measurement  errors  are  transferred  from  a  known 
calibration  target  to  the  target  under  test.  Measurements  are  made  of  the  fields  scattered 
by  the  target  and  its  mount,  E |>,  the  target  mount  alone,  M ,  a  calibration  object  and 
its  mount,  E £,  and  the  calibration  object  mount  alone,  EqM.  A  calibrated  measurement 
of  the  target,  aj  is  obtained  with  vector  background  subtraction  using  the  formula, 


Gyp  = 


ET  -  etm 


2 


&cal  i 


(i.i) 


where  <rcai  is  the  calculated  RCS  of  the  simple  calibration  target  used.  The  quantity 
Ej>  —  Ejpm  ls  assumed  to  be  a  clean  measurement  of  the  target  alone.  Likewise,  the 
quantity  Eq  —  EqM  is  assumed  to  be  a  clean  measurement  of  the  calibration  object  alone. 
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A  perfectly  conducting  sphere  with  its  Mie  series  solution  is  commonly  used  as  a 
calibration  object.  Other  perfectly  conducting  bodies  of  revolution  (BOR)  may  also  be 
used — and  may  produce  a  better  calibration.  For  example,  the  geometry  of  a  right-circular 
cylinder  eliminates  certain  interactions  between  the  itself  and  the  pylon  on  which  it  is 
mounted.  Use  of  other  BORs,  however,  is  dependent  upon  the  ability  to  calculate,  or 
otherwise  determine,  the  expected  scattering  to  within  a  specified  tolerance. 

Current  computer  codes  which  calculate  the  scattering  from  bodies  of  revolution 
implement  the  Method  of  Moments.  They  subdivide  the  scattering  surface  into  sections 
and  use  piecewise  expansions  of  the  unknown  current  density  in  an  integral  equation  that 
describes  the  scattering  while  satisfying  boundary  conditions.  The  method  of  expanding 
the  current  on  sub-sections  of  the  surface  results  in  two  deficiencies.  First,  a  prohibitively 
large  number  of  these  expansion  functions  are  required  to  accurately  describe  the  current 
density  on  electrically  large  bodies.  Second,  if  sufficient  accuracy  is  not  obtained,  a  re- 
sectionalizing  and  re-calculation  of  the  current  density  on  the  surface  is  required,  resulting 
in  lost  computational  time. 

The  approach  proposed  here  is  to  expand  the  unknown  current  in  the  equations  using 
functions  which  have  the  entire  length  of  the  surface  as  their  domain.  Functions  are  chosen 
that  better  represent  the  expected  current  density  on  the  surface.  The  continuous  nature 
of  entire-domain  functions  allows  them  to  better  describe  the  current  density  on  a  smooth 
surface  such  as  that  of  a  BOR.  Also,  since  re-discretization  is  not  required  to  add  additional 
expansion  terms,  accuracy  can  be  improved  without  discarding  previous  calculations.  Both 
of  the  shortcomings  described  for  subdomain  sectionalization  are  avoided. 

The  goal  of  this  research  is  to  demonstrate  the  utility  of  entire-domain  basis  func¬ 
tions  in  computing  the  scattering  from  bodies  of  revolution.  Chapter  II  briefly  outlines  the 
Moment  Method  which  is  used  to  solve  the  radar  scattering  equations.  It  also  describes 
subdomain  basis  functions  that  are  commonly  used  to  expand  the  unknown  current  density 
in  the  Moment  Method  solution.  Entire-domain  basis  functions  are  introduced  and  the 
Chebyshev  polynomials  are  chosen  as  a  suitable  family  of  expansion  functions.  Advan¬ 
tages  and  disadvantages  are  discussed  for  each  type  of  function.  An  electric-field  integral 
equation  is  developed  in  Chapter  III  for  scattering  from  an  arbitrary  body  of  revolution. 
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The  Chebyshev  polynomials  are  then  used  in  Chapter  IV  as  the  basis  for  expanding  the 
unknown  current  density  induced  on  a  perfectly  conducting  sphere  by  plane  wave  exci¬ 
tation.  The  resulting  matrix  equations  are  implemented  in  a  computer  code  and  solved. 
The  data  reported  in  Chapter  V  demonstrates  that  the  RCS  can  be  computed  in  fewer 
entire-domain  basis  expansion  modes  than  subdomain  basis  functions  under  the  same  con¬ 
ditions.  An  analysis  of  solution  matrix  fill  time  and  potential  savings  is  also  performed. 
Chapter  VI  summarizes  the  results.  Although  it  is  shown  that  fewer  entire-domain  modes 
are  required  for  the  RCS  computation,  the  singular  nature  of  the  integrals  adds  to  the 
complexity  of  their  robust  implementation.  The  method  of  equivalent  separation  of  the 
testing  functions  and  basis  functions  used  in  subdomain  methods  did  not  yield  an  accurate 
solution  for  spheres  of  small  electrical  size.  A  more  rigorous  handling  of  the  singularity 
and  other  recommendations  for  further  research  are  suggested.  The  computer  code  and 
additional  data  are  included  in  Appendix  A  and  B,  respectively. 
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II.  Background 

The  problem  of  improving  the  accuracy  and  computation  time  of  calculating  the  scattering 
from  bodies  of  revolution  centers  on  the  numerical  solution.  An  important  step  in  a 
numerical  solution  is  the  choice  of  an  appropriate  basis  function  set  by  which  to  expand 
the  unknown  quantity.  This  chapter  discusses  the  numerical  method  that  will  be  used  to 
solve  the  integral  equations  that  describe  the  scattering  from  a  body  of  revolution  caused  by 
an  incident  plane- wave.  It  also  outlines  in  some  detail  the  various  types  of  basis  functions 
and  reasons  for  their  consideration.  Finally,  justification  for  the  choice  of  entire-domain 
basis  functions  is  given. 

2.1  Numerical  Solution  by  the  Method  of  Moments 

While  several  methods  exist  for  determining  an  asymptotic,  or  high  frequency,  ap¬ 
proximate  solution  to  a  scattering  problem,  the  task  here  is  to  find  a  solution  that  converges 
to  the  exact  solution  in  regions  where  the  asymptotic  methods  fail.  One  such  method,  the 
Method  of  Moments,  produces  a  solution  that  can  converge  to  the  exact  solution.  This 
method  is  the  foundation  for  many  variants  of  differing  names. 

The  Method  of  Moments  was  first  developed  to  solve  elastodynamic  problems  and 
later  generalized  to  a  method  of  solution  for  a  wide  variety  of  problems  involving  linear 
integral  and  differential  operations  [14].  Also  called  the  Moment  Method,  it  is  a  numer¬ 
ical  technique  by  which  a  linear  operator  equation — one  in  which  the  unknown  quantity 
is  imbedded  in  an  integrand  or  the  argument  of  a  differentiation  and  cannot  be  solved 
for  analytically — can  be  solved  approximately  by  transforming  the  operator  equation  into 
a  system  of  linear  equations.  Once  the  problem  is  formulated  by  this  method,  the  solu¬ 
tion  comes  by  computing  the  inverse  matrix  and  performing  a  matrix  multiplication.  An 
integro-differential  equation  such  as  the  one  that  will  be  developed  for  plane-wave  scatter¬ 
ing  from  bodies  of  revolution  can  be  written  as  a  linear  operator  acting  upon  the  unknown 
current  density  to  produce  a  functional  effect.  Such  an  operation  is  written  mathematically 
as 

/  =  ![«],  (2.1) 
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where  L  is  the  linear  operator,  u  is  the  unknown  function,  and  /  is  the  effect  which  is  often 
called  the  forcing  function.  The  unknown  quantity  is  expanded  into  a  series  of  terms.  Each 
term  is  one  of  a  set  of  linearly  independent  functions — called  basis  functions — multiplied 
by  a  coefficient  which  describes  that  term’s  relative  importance  to  the  total  sum.  However, 
even  with  computers,  an  infinite  number  of  terms  cannot  be  handled  computationally,  so 
the  series  is  truncated  to  a  finite  series  of  terms  plus  a  final  term  which  represents  all 
remaining  terms  (the  error,  €), 

oo 

«  =  £  Infin  (2.2) 

n=— oo 
+W2 

=  E  Wn  +  e.  (2.3) 
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Or,  writing  as  a  one-sided  sum, 
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Substituting  the  above  expression  for  u  into  Equation  (2.1),  the  constants  can  pass  through 
the  operator  because  of  its  linearity  giving 
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where  re  is  called  the  residual.  The  residual  is  then  minimized  by  taking  an  inner  product 
with  a  family  of  testing  functions,  0m,  and  forcing  it  to  be  orthogonal  to  the  testing 
function.  Galerkin  does  this  by  choosing  V’n  =  [11].  For  many  functions  involving 

complex  calculations,  a  suitable  inner  product  is 


(fim  i  tyn)  —  J  j  @m  '  ‘finds. 
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Applying  this  inner  product  to  Equation  (2.6)  and  choosing  to  enforce  (0m,re)  —  0  results 
in  a  system  of  linear  equations, 

N 

D  In  ( L[9m ,  i’n])  ~  (9m,  />  =  <>,  m  =  0,  1, 2, . . . ,  N.  (2.8) 

n=0 

The  choice  of  N  testing  functions  and  N  basis  functions  leads  to  a  system  of  N  equations 
with  N  unknowns,  In.  Rearranging  gives 

N 

<W>=  'E'WemMh  m  =  0)  1,2, (2.9) 
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which  can  be  easily  written  in  simple,  block  matrix  notation 

IVIjvxi  =  [Z]nxN  [-Hjvxi.  (2.10) 


where 


Win  =  (°mJ) 

W]mn  =  ( I  \@m  i  V’n]) 

[I]  n  =  ^ 

for  m=  0,  1,  2,  . . N 

n=  0,  1,  2,  ....  N.  (2.11) 

When  solved  for  the  unknown, 

WiVxi  =  [Z&'xiV  m,x„  (2-12) 

where  the  elements  of  [/]  are  the  coefficients  of  the  expansion  terms  of  the  unknown  quantity 
in  the  original  operator  equation.  The  accuracy  of  the  numerical  expansion  is  limited  by 
the  number  of  expansion  terms,  N. 
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2.2  Choosing  a  Basis  Function 

Although  there  are  infinitely  many  functions  which  could  be  used  as  basis  func¬ 
tions,  an  illogical  choice  of  functions  may  limit  the  accuracy  that  can  be  obtained  from 
the  solution  or  may  demand  a  prohibitively  large  number  of  functions  for  convergence. 
An  appropriate  basis  set  for  calculating  the  scattered  electromagnetic  fields  from  a  body 
contains  functions  that  approximate  the  suspected  current  density  on  the  surface  of  the 
body  (rapid  convergence),  and  have  inter-function  relationships  that  offer  numerical  ad¬ 
vantages  in  the  computation  process  (saved  computation  time).  It  should  be  noted  that 
some  functions  may  meet  the  requirements  for  a  basis  set,  but  may  never  be  able  to  ap¬ 
proach  the  exact  solution  with  accuracy.  This  situation  would  occur,  for  example,  if  one 
were  to  attempt  a  solution  by  using  basis  functions  that  have  smoother  properties  than  the 
unknown  being  represented  [1] .  Factors  affecting  the  choice  of  a  basis  function  set  involve 
the  desired  solution  accuracy,  the  relative  complexity  of  the  resulting  matrix  entries,  and 
the  computational  constraints  (resources  and  time)  that  limit  the  ultimate  matrix  size  [11]. 

2.2.1  Subdomain  Basis  Functions.  One  of  the  simplest,  and  therefore  most 
common,  basis  sets  that  are  used  are  sets  of  subdomain  functions.  That  is  to  say  they  are 
functions  that  are  defined  only  over  the  domain  of  consecutive  sections  of  the  structure. 
A  surface  is  discretized  by  dividing  it  into  N  sections.  The  sections  need  not  be  of  equal 
size  but  commonly  are  for  computational  simplicity.  A  function  is  defined  over  the  limits 
of  each  segment.  Simple  functions  may  be  a  pulse  (constant  amplitude),  a  triangle,  an 
arc  of  a  sinusoid,  etc.  The  unknown  current  or  other  quantity  being  calculated  is  then 
approximated  by  multiplying  each  piecewise  function  by  an  amplitude  factor.  This  is 
described  mathematically  by  the  Equation  (2.13)  where  f(x)  is  the  unknown  function 
being  approximated,  n  is  an  index  for  each  section,  an  is  the  amplitude  factor  and  gn(xf) 
is  the  basis  function, 

f(x)  ~^r,an9n(x').  (2.13) 
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Given  the  subdivision  of  an  axis  shown  in  Figure  2.1,  pulse-or  piecewise  constant-basis 
functions  are  placed  on  each  section.  The  pulse  function  is  defined  by 


9n{x') 


r 

1  %n— 1  —  ®  —  xn 

< 

0  elsewhere 


(2.14) 


The  height  of  each  pulse  corresponds  to  the  an  for  that  section.  The  curve,  /(#),  repre¬ 
sented  by  the  semi-circle  is  approximated  by  the  pulses.  An  example  is  Glisson’s  use  of 
pulse  basis  functions  to  expand  the  current  in  his  body  of  revolution  solution  [4]. 

Triangular — or  piecewise  linear — basis  functions  are  similarly  used,  but  they  are  usu¬ 
ally  defined  with  a  domain  that  covers  two  consecutive  segments.  They  are  used  in  an 
overlapping  fashion  as  shown  in  Figure  2.2  so  that,  when  summed,  the  amplitude  approx¬ 
imates  the  desired  curve.  Triangular  expansion  functions  are  useful  in  that  they  generally 
simplify  the  evaluation  of  the  integro-differential  operator  by  approximating  the  curve  over 
each  segment  as  linear.  Both  Mautz’  formulation  [8]  and  Rogers’  JRMBOR  [13]  are  exam¬ 
ples  of  the  use  of  triangular  basis  functions.  The  triangle  basis  functions  approximate  the 
unknown  curve  with  straight  lines  whose  derivatives  are  constants.  Many  moment  method 
based  computer  codes  use  triangular  basis  functions  because  they  produce  a  continuous 
approximation  to  the  current. 

Subdomain  expansion  functions  are  not  advantageous  when  it  comes  to  adding  more 
terms  to  the  expansion  to  increase  accuracy.  Prior  results  based  upon  an  N-segmented, 
equal-interval  domain  must  be  discarded  to  re-discretize  the  domain  into  greater  than  N 
segments. 


2.2.2  Entire  Domain  Basis  Functions .  Basis  functions  may  also  be  defined  over 
the  entire  length  of  a  surface  under  consideration.  Such  functions  are  appropriately  called 
entire  domain  basis  functions.  Equation  (2.13)  is  still  valid  with  an  appropriate  entire- 
domain  definition  of  gn(x ').  Since  the  structure  is  not  segmented  as  in  the  case  of  sub- 
domain  basis  functions,  each  expansion  term,  n,  is  referred  to  as  a  mode  rather  than  a 
section  or  a  segment. 
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Figure  2.2  Triangular  (piecewise  linear)  subdomain  basis  functions. 
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A  familiar  example  of  entire  domain  basis  functions  is  the  Fourier  series  in  which  a 
function  is  expanded  in  terms  of  exponentials  (or  sines  and  cosines) , 
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(2.15) 

cn  =  r  x(<t>)e-jn*d<f>. 
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(2.16) 

An  electrical  pulse  can  be  approximated  in  this  way  by  adding  together  higher  and  higher 
orders  of  sinusoids  from  an  oscillator  until  the  edges  of  the  pulse  are  sufficiently  steep  and 
the  center  ripple  is  sufficiently  small  for  a  given  application.  Balanis  [2]  points  out  that  an 
entire  domain  sine  is  an  appropriate  basis  for  the  current  on  a  dipole  antenna  because  the 
current  is  known  to  be  zero  at  the  ends  and  to  oscillate  approximately  as  a  sine  between 
the  ends  according  to  the  frequency  of  the  excitation.  Poulsen  [12]  has  experimented  with 
sine  and  cosine  basis  functions  for  computing  the  fields  induced  on  a  frequency-selective 
surface  of  crossed-dipole  arrays.  While  no  evidence  of  any  other  use  of  entire-domain  basis 
functions  for  computing  the  scattering  from  bodies  of  revolution  was  found,  Kitazawa  [7] 
and  Butler  [3]  used  Chebyshev  polynomials  to  solve  for  the  fields  resonant  in  illuminated 
rectangular  slots  and  waveguides. 

Unlike  subdomain  basis  functions,  which  when  the  interval  is  divided  into  sufficiently 
many  intervals  can  approximate  an  arbitrary  curve  with  the  appropriate  amplitudes  ap¬ 
plied,  entire-domain  modes  can  only  be  used  effectively  when  the  general  nature  of  the 
curve  is  known  a  priori  allowing  choice  of  a  basis  that  exhibits  similar  features  [1].  This 
can  be  a  disadvantage  or  an  advantage  depending  upon  the  extent  of  information  known  or 
assumed  about  the  problem  at  hand.  A  clear  advantage  of  entire-domain  basis  functions 
is  that  if  accuracy  is  not  sufficient  for  an  initial  finite  number  of  modes,  additional  modes 
can  be  added  to  the  expansion  while  retaining  the  sum  of  previous  modes. 

2.2.3  Physical  Basis  Functions.  One  subset  of  entire  domain  basis  functions  is 
physical  basis  functions.  These  functions  are  unique  in  that  they  are  based  upon  knowledge 
of  current  components  that  have  actual  physical  interpretations.  For  example,  some  bodies 
exhibit  scattering  currents  that  can  be  traced  to  actual  physical  features  or  properties  of 
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the  body.  Kent  [6]  found  them  useful  in  reducing  matrix  size  in  the  solution  of  scattering 
from  a  dielectric  coated  strip.  Physical  basis  functions  are  particularly  useful  for  bodies  of 
larger  electrical  size  because  they  can  reduce  the  order  of  the  system  to  be  solved. 

2.2.4  Hybrid  Basis  Functions .  Some  methods  may  use  several  different  types  of 
basis  functions  simultaneously.  One  hybrid  method  involves  using  basis  functions  which 
represent  the  current  as  determined  by  a  high  frequency  asymptotic  approach  such  as 
Physical  Optics  or  Geometrical  Theory  of  Diffraction  and  combining  it  with  subdomain 
basis  functions  near  edges  or  other  singularities  [15].  This  approach  also  has  as  its  goal 
reducing  the  order  of  the  system  of  equations  to  be  solved  which  has  benefits  for  solutions 
to  problems  involving  electrically  large  bodies. 

Combining  physical  basis  functions  with  other  basis  functions  may  also  be  useful.  A 
large  body  may  be  analyzed  with  physical  basis  functions  and  subdomain  basis  functions, 
for  example.  Physical  basis  functions  are  used  to  reduce  the  number  of  unknowns  over 
electrically  large  regions,  and  subdomain  basis  functions  are  used  to  increase  accuracy 
near  edges  or  other  smaller  features. 

2.2.5  Vector  Basis  Functions.  The  subdomain  and  entire-domain  basis  functions 
as  described  above  have  only  employed  a  scalar  representation  which  can  be  used  in  the  so¬ 
lutions  to  three-dimensional  problems  by  separately  expanding  the  coordinate  components 
of  the  unknown  function.  Under  some  conditions,  this  may  not  be  desirable.  If  a  problem 
centers  on  the  solution  at  interfaces  between  different  media,  a  boundary  condition  can 
be  difficult  to  impose  on  an  expansion  that  does  not  naturally  separate  into  components 
along  the  boundary.  Also,  scalar  expansions  tend  to  be  continuous,  while  with  some  vector 
quantities,  enforcing  this  characteristic  can  introduce  additional  errors.  Peterson  shows 
that  vector  basis  functions  provide  not  only  a  function  amplitude  in  a  region,  but  a  di¬ 
rection  as  well  [11].  Figure  2.3  shows  a  notional  example  of  directional  functions  within 
a  triangular,  multi-dimensional,  subdomain  cell.  The  vector  basis  functions  can  be  de¬ 
fined  appropriately  at  medium  interfaces  to  satisfy  continuous  or  discontinuous  boundary 
conditions. 


2-9 


2.3  Testing  Functions 

In  order  to  determine  the  electromagnetic  scattering  from  a  body  of  revolution,  the 
body  is  illuminated  by  an  incoming  wave  and  measurements  are  made  at  various  points 
around  the  body  to  determine  the  scattering  pattern.  Accuracy  depends  upon  the  number 
of  measurements  made  as  well  as  the  type  of  measuring  device  used.  This  measurement  is 
accomplished  numerically  by  means  of  testing  functions  placed  around  the  body.  Just  as  in 
the  case  of  the  basis  functions  used  to  expand  the  unknown  current,  the  testing  functions 
can  be  subdomain  or  entire-domain  functions.  The  measurement  of  the  scattering  by  each 
testing  function  mode  can  be  thought  of  as  a  passive  measurement  of  the  scattering  of  the 
basis  function  modes  which  is  induced  by  the  incident  wave.  Electromagnetic  reciprocity 
also  allows  thinking  about  this  as  the  same  problem  in  reverse.  In  this  case,  the  body 
is  excited  by  each  test  function  which  results  in  scattering  by  each  basis  function.  The 
resulting  fields  are  then  measured.  It  is  from  this  point  of  view  that  the  numerical  solution 
is  developed.  The  reaction  of  the  basis  functions  to  each  testing  function  is  computed  via 
an  inner  product. 

Any  number  of  kinds  of  testing  functions  exist.  However,  two  kinds  are  commonly 
used  because  of  the  numerical  advantages  they  offer.  The  first,  and  possibly  the  most 
commonly  used,  type  is  the  Dirac  delta  distribution.  In  an  integral  equation,  these  func¬ 
tions  allow  elimination  of  one  dimension  of  integration  for  each  dimension  in  which  they 
are  used,  thereby  simplifying  the  numerical  computation.  For  this  simplicity,  a  trade-off 
is  made  by  relaxing  the  boundary  conditions  so  that  they  are  only  enforced  at  discrete 
points.  Therefore,  proper  placement  of  the  points  is  important. 

Another  possible  choice  of  testing  functions  is  a  set  of  functions  identical  to  the  expan¬ 
sion  basis  functions.  This  is  known  as  Galerkin’s  method.  For  certain  integral  equations, 
choosing  testing  functions  equal  to  the  basis  functions  allows  an  orthogonal  projection  of 
the  unknown  quantity  onto  the  basis  functions  resulting  in  numerical  convergence  in  the 
limit  as  the  number  of  basis  functions  goes  to  infinity.  The  integral  operators  that  arise 
in  electromagnetics  are  generally  not  quite  so  ideal.  However,  Galerkin’s  method  is  of¬ 
ten  chosen  for  other  reasons,  one  being  that  its  use  with  electrical  field  integral  equations 
produces  matrices  with  diagonal  symmetry.  This  permits  an  approximately  50%  reduc- 
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tion  in  matrix  fill  computation  time.  This  choice  is  acceptable  in  electromagnetic  integral 
equations  because  even  in  the  case  where  orthogonality  is  possible,  an  infinite  number  of 
functions  can  never  be  used.  An  error  always  exists  and  is  monitored  and  controlled  by 
the  number  and  type  of  functions  used  [11]. 

2.4  Summary  and  Implications  to  Bodies  of  Revolution 

The  types  of  basis  functions  that  exist  are  as  varied  as  the  types  of  problems  to 
be  solved.  The  categories  include  subdomain  functions,  entire-domain  functions,  hybrid 
basis  functions,  and  directional  (vector)  functions.  The  problem  of  interest  here  is  one  of 
computing  the  scattering  from  a  body  of  revolution  in  free  space.  While  subdomain  basis 
functions  are  commonly  used,  they  are  not  the  most  time-efficient  choice  for  a  desired 
accuracy.  In  this  development,  the  surface  shape  of  the  body  is  kept  arbitrary  throughout 
the  solution,  eliminating  the  usefulness  of  physical  basis  functions.  Also,  since  this  problem 
is  posed  as  a  free-space  problem  (constant  medium) ,  vector  basis  functions  are  not  necessary 
or  helpful.  Choosing  entire-domain  basis  functions  offers  several  desired  advantages.  Using 
entire-domain  basis  functions  results  in  a  smoother  approximation  to  the  actual  current 
density  than  using  subdomain  basis  functions.  An  accurate  approximation  of  the  current 
density  is  also  possible  in  fewer  modes  assuming  an  appropriate  entire-domain  function  is 
used  (exact  convergence  is  only  possible  if  the  basis  function  is  also  an  exact  solution). 
Lastly,  computational  time  is  saved  due  to  the  ability  to  incrementally  monitor  and  add 
modes  as  required  without  re-computation  of  previous  modes.  In  using  entire-domain  basis 
functions,  proper  placement  of  delta  testing  functions  is  not  obvious  and  Galerkin’s  method 
is  used.  The  integral  equation  derived  in  Chapter  III  is  written  as  a  matrix  equation  in 
Chapter  IV  for  an  arbitrary  choice  of  basis  and  testing  functions.  Chebyshev  entire-domain 
basis  functions  and  Galerkin  testing  are  used  to  generate  the  results  in  Chapter  V. 
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Ill .  Theory 

In  this  chapter,  the  geometry  for  an  arbitrary  perfectly  electric  conducting  (PEC)  body 
of  revolution  is  established.  For  a  PEC  body,  the  electric  field  componenets  tangential  to 
the  surface  must  vanish  at  every  point  on  the  surface.  Using  this  boundary  condition,  an 
electric  field  integral  equation  is  formed  to  describe  the  scattering  from  the  body.  Also, 
the  Chebyshev  polynomials  that  will  be  used  as  basis  functions  in  the  numerical  solution 
of  the  integral  equation  are  described.  The  conventions,  notation,  and  overall  development 
follows  that  of  Mautz  and  Harrington  [9]. 

3.1  Definition  and  Geometry 

A  body  of  revolution  is  created  by  rotating  a  generating  arc  around  an  axis.  Consider 
a  generating  arc  defined  by  a  continuous  function  in  the  XZ- plane  whose  endpoints  are 
on  the  £-axis.  The  function  may  be  defined  in  piecewise  fashion,  and  does  not  need  to 
be  differentiable  in  the  sense  that  the  resulting  curve  can  have  corners.  This  function  is 
then  rotated  axially  around  the  z-axis  to  create  a  closed  surface  that  is  symmetric  in  the 
cylindrical  direction  <j>.  The  surface  and  the  volume  it  encloses  is  called  a  body  of  revolution 
because  it  is  created  by  revolving  a  function  around  an  axis,  not  because  of  any  motion 
of  the  body  itself.  Such  a  body  of  revolution  is  depicted  in  Figure  3.1.  Although  the 
body  is  depicted  as  being  defined  entirely  in  the  positive-z  half-space,  this  is  not  necessary. 
However,  the  equations  used  herein  to  define  the  geometric  parameters  of  the  body  are 
written  assuming  positive-z  and  may  not  be  valid  for  other  conventions. 

It  is  convenient  to  write  the  generating  arc  as  a  parametric  equation  of  a  curvilinear 
coordinate  system  on  the  surface  of  the  body  of  revolution.  In  Figure  3.1,  p,  and  z  are 
the  ordinary  cylindrical  coordinates  and  i  and  (f>  are  orthogonal  unit  vectors  at  a  point, 
5,  on  the  surface  which  satisfy  <j>  X  t  =  h.  Letting  to  be  the  zero  arclength  point  of  the 
generating  arc  (located  at  the  origin),  i  points  in  the  direction  of  increasing  arclength, 
and  <j>  points  in  the  direction  of  increasing  cylindrical  <f>.  This  leads  to  choosing  n  as 
the  outward  normal  unit  vector  to  the  surface  S.  The  body  of  revolution  is  excited  by 
an  incident  plane  wave  as  shown  in  Figure  3.2.  The  vector  kf  indicates  the  direction  of 
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Figure  3.1  Geometry  of  the  body  of  revolution 


3-2 


plane  wave  propagation.  The  coordinates  <j)t  and  Of  indicate  the  transmitter  location.  The 
vectors  <j>t  and  0t  are  unit  vectors  in  the  </>t  and  6t  directions,  respectively.  The  resulting 
scattered  field  components  will  be  measured  in  the  monostatic  (backscatter)  direction.  It 
will  be  assumed  for  the  purposes  of  this  problem  that  the  body  of  revolution  has  a  perfect 
electrically  conducting  surface. 

3.2  Integral  Equations 

The  problem  of  scattering  from  the  body  of  revolution  is  solved  as  an  electromagnetic 
boundary  value  problem  by  forming  either  an  electric  field  integral  equation  or  a  magnetic 
field  integral  equation.  It  has  been  shown  by  Mautz  and  Harrington  [9]  that  the  solutions 
to  the  E-field  and  H-field  integral  equations  deteriorate  near  internal  resonances  of  the 
conducting  surface,  but  a  solution  involving  a  weighted  linear  combination  of  the  E-field 
and  H-field  integral  equations  does  not.  However,  they  also  show  that  although  the  error  in 
current  density  near  the  resonances  using  the  E-field  integral  equations  can  be  tremendous, 
the  error  in  scattering  cross-section  is  quite  small  compared  to  that  of  the  H-field  solution 
and  is  comparable  to  that  of  the  combined  field  solution.  Using  both  E-field  and  H-field 
integral  equations  in  a  combined  equation  effectively  doubles  the  computation  time  for 
each  element  in  the  solution  while  giving  minimal  improvement  in  the  scattering  cross- 
section  data.  For  this  reason,  development  and  use  of  the  E-field  integral  equation  for 
this  concept  demonstration  is  considered  sufficient.  The  H-field  integral  equation  can  be 
similarly  derived  and  used  to  trade  computation  time  for  increased  accuracy. 

Maxwell’s  equations  require  that  the  tangential  electrical  fields  at  the  surface  of  a 
perfect  electrical  conductor  sum  to  zero.  This  boundary  condition  is  written  as 

C  +  C  =  0,  (3.1) 

where  Emc  is  the  incident  electric  field  and  Es  is  the  field  scattered  from  the  surface.  The 
subscript  tan  denotes  the  tangential  component  on  S .  The  incident  field  is  defined  to  be 
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either  9  ox  <j>  directed  and  is  written  as 


Elnc  =  9kr)e-jkt-f 

(3.2) 

Efc  =  ^krje-jk‘-f, 

(3.3) 

where  k  is  the  propogation  constant,  r/  is  the  intrinsic  impedance  (377Q  for  free-space), 
and  f  is  the  radius  vector  from  the  origin  to  a  point  on  the  surface.  The  scattered  field  is 
found  via  the  magnetic  vector  potential,  A  and  the  electric  scalar  potential,  $e, 


ES  =  -jojA(J)  -  V#e(J). 


(3.4) 


The  magnetic  vector  potential,  A ,  and  the  electric  scalar  potential,  4>e,  written  in  terms 
of  the  free-space  Green’s  function  are  [2] 


rjds' 

(3.5) 

,jk\r-r’\ 

- — —\ds’. 
r  —  r 

(3.6) 

Here,  f  and  f'  are  vectors  to  the  measurement  and  source  points,  respectively.  The  quantity 
J(f')  is  the  electric  current  density  on  5,  k  is  the  propagation  constant,  and  p  and  e  are 
the  permeability  and  permittivity,  respectively.  The  operator  V-  is  the  divergence.  The 
quantity  |  r  —  r '  |  can  be  written  in  terms  of  the  cylindrical  coordinates  as 


R=\ 


-/ 

r 


|=  \  ( p  -  pf )2  +  (z-  z')2  +  4/y/sin1 


(3.7) 


where  z ,  and  v  are  functions  of  arclength,  t  as  shown  in  Figure  3.1.  Similarly,  p\  z', 
and  vf  are  functions  of  tf.  Substituting  Equation  (3.4)  with  (3.5),  (3.6)  and  (3.7)  into 
Equation  (3.1)  gives 


TPinc 
& tan 


flsm 


pjkR 


R 


-ds'  + 


juiAne 


Jjs[v.J(P)\ e 


,jkR 


R 


-ds'  =  0, 


(3.8) 
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which  is  an  electric  field  integral  equation  describing  the  tangential  electric  fields  at  the 
surface.  This  is  a  complicated  integro-differential  equation  where  J  is  the  unknown  quantity 
which  can  be  determined  numerically  using  the  Method  of  Moments.  Once  J  is  known,  it 
can  be  used  to  compute  the  scattered  electric  fields. 

3.3  Chebyshev  Polynomial  Basis  Functions 

There  are  many  entire-domain  functions  that  could  be  used  as  a  basis  to  expand  the 
unknown  current  density  on  the  surface  of  the  body.  Some  common  orthogonal  functions 
to  which  the  entire-domain  of  the  generating  arc  could  be  mapped  are  sinusoids,  Jacobi 
polynomials,  Chebyshev  polynomials,  Legendre  polynomials,  Maclaurin  series,  Laguerre 
polynomials,  and  Hermite  polynomials  [10]. 

Since  one  of  the  purposes  in  using  an  entire-domain  basis  is  to  reduce  the  number 
of  modes  or  terms  required  to  approximate  the  current  density  when  expanded  as  a  se¬ 
ries,  the  shape  of  the  basis  functions  themselves  should  closely  model  the  functionality  of 
the  expected  current  density.  The  function  must  also  have  higher  modes  that  exhibit  a 
sufficient  degree  of  oscillatory  behavior  or  other  fluctuation  in  its  derivatives  to  define  the 
changes  in  the  current  density  across  the  length  of  the  body.  The  current  density  should 
be  systematically  improved  by  increasing  the  number  of  basis  functions. 

The  Chebyshev  polynomials  of  the  first  kind  are  found  to  have  the  required  char¬ 
acteristics  in  the  domain  of  —  1  <  s  <  1.  An  additional  advantage  of  this  choice  is  that 
the  positive  and  negative  ends  are  symmetric  and  finite.  Table  3.1  gives  the  functional 
form  of  the  first  nine  polynomials  which  are  illustrated  in  Figure  3.3.  The  polynomials 
are  easily  referenced  in  a  computer  code  by  a  lookup  table  of  the  coefficients.  The  higher 
order  polynomials  can  be  generated  by  a  recursion  relationship  [10], 


Tn{x)  —  2,xTn—\(x)  Tn_2(x). 


(3.9) 
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arclength 


arclength 


arclength 


Figure  3.3  The  first  nine  modes  of  the  Chebyshev  polynomial  of  the  first  kind  (T0-T8) 


To  =  1 
T\  =  x 
T2  =  2a;2 -1 
To  =  4a;3  —  3a; 

T4  =  8a;4  —  8a;2  +  1 
T5  =  16a;5  -  20a;3  +  5a: 

T6  =  32a;6  -  48a;4  +  18a;2  -  1 
T7  =  64a;7-  112a;5 +  56x3-  lx 
T8  =  128x8  -  256a;6  +  160a:4  -  32a;2  +  1 

Table  3.1  Chebyshev  polynomials  of  the  first  kind. 


The  Chebyshev  basis  functions  are  orthogonal  with  respect  to  a  weighting  function,  w(x)  — 

1/(1  -*2)1/2, 


j\m(x)Tn(x)( l-x2)~1/2  dx  = 


0, 

71-/2, 

*■> 


ra  /ft 

m  =  n^0  • 

ra  =  n  =  0 


(3.10) 


However,  the  weighting  is  not  necessary  for  convergence  if  the  operator,  L,  is  self-adjoint 
with  respect  to  the  inner  product,  (La,  6)  =  (a,  L6),  and  is  positive  definite,  (a,  La)  >  0 
for  all  nonzero  a  [11].  The  linear  independence  of  the  Chebyshev  polynomials  is  sufficient 
to  force  the  residual  error  of  the  inner  product  to  zero  and  establishes  the  convergence  of 
the  finite  expansion.  In  practice,  Galerkin  testing  is  used  with  the  electric-field  integral 
equation  without  regard  to  true  orthogonality  because  it  produces  matrices  with  diagonal 
symmetry  reducing  the  required  computations. 

In  summary,  an  electric-field  integral  equation,  Equation  (3.8)  was  formed  as  a  func¬ 
tion  of  the  unknown  current  density.  The  equation  applies  on  the  surface  and  is  based  on 
the  vanishing  tangential  electric  fields.  Also,  the  Chebyshev  polynomials  were  defined  and 
will  be  used  to  expand  the  current  density  in  the  numerical  solution  outlined  in  Chapter  IV . 
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IV.  Application  to  the  Integral  Equation 

In  this  chapter,  the  electric  field  integral  equation  given  in  Equation  (3.8)  is  solved  for  the 
unknown  current  density  by  means  of  the  Method  of  Moments  using  entire-domain  basis 
functions. 

Equations  (3.1)  and  (3.4)  are  rewritten  in  terms  of  the  tangential  component  of  the 
incident  field, 


Tpmc  _ 

^tan  *— 


-Ef 


tan 


=  -[-V$e(J)-juA(J)] 


tan  ' 


Define  a  linear  operation  on  the  unknown  current  density,  J,  as 


L(J)  =  \juA(J)  +  V*(J)]taw. 


(4.1) 


(4.2) 


The  current  density  is  written  in  a  series  expansion  of  linearly  independent  basis  functions 

j  =  (4-3) 

3 

Because  the  unknown  currents  J  are  periodic  in  <f>  (body  of  revolution),  each  ipj  can  be 
represented  as  a  Fourier  series  in  <f>.  The  Fourier  series  is  an  entire-domain  expansion  where 
the  domain  is  0  <  <f>  <  2n  and  the  basis  is  the  complex  exponential  ejn<?5>  which  represents 
a  Fourier  mode  indexed  by  n.  Now, 

J{t'A ')  (4.4) 

n  j 

Primes  are  used  here  and  in  the  remainder  of  this  document  to  denote  association  with 
a  source  (surface  scattering)  point.  Unprimed  quantities  refer  to  field  points  where  the 
electric  field  is  measured  or  tested.  The  current  density  at  each  point,  t,  along  the  gen¬ 
erating  arc  is  completely  described  by  the  sum  of  the  current  densities  in  two  orthogonal, 
curvilinear  component  directions,  t  and  O.  Separately  expanding  InjJj  in  each  direction 
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gives 


Inj  Jj 


01 + 44 
Km +Km, 


(4.5) 


where  t'  is  a  point  on  the  generating  arc,  0  <  t*  <  tjy  and  fj(t')  is  the  j’th  basis  function. 
The  current  expansions  are,  therefore,  defined  as 


4  = 

4  =  vmm.  (4.6) 


Substituting  Equation  (4.4)  and  (4.5)  into  Equation  (4.3)  gives  a  final  expression  for  the 
current  density, 

j = EE  [Km + Km]  m.  (4.7) 

n  j 

The  integral  equations  now  can  be  written  as 


Elan  =  L 


Y^  p*  p* ,  _l  Y^  T&'  j^[ 

/  v  Ti3unj  1  /  _j  njunj 


n,3 


n,3 


(4.8) 


The  linearity  of  L  allows  interchanging  the  order  of  summation  and  integration/differentiation, 
for  finite  sums,  to  produce 


4  =  E  Kl  (K)  +  41  (4)]  ■  (4-9> 

n,j 

Since  each  Fourier  expansion  mode  is  linearly  independent  from  and  orthogonal  to  every 
other,  the  above  formulation  is  separable  in  n  and  involves  a  separate  solution  to  determine 
the  current  coefficients  for  each  Fourier  mode.  Rewriting  as  a  separate  equation  for  each 
n  permits  simplification  of  the  notation  by  dropping  the  subscript  n, 

e=E  lIi'L{Ji)+IfL  {■>!')]•  v » •  (44°) 
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Testing  functions  are  chosen  equivalent  to  the  basis  (current  expansion)  functions 
according  to  Galerkin’s  method, 


WL  =  *m+ 

W*i  =  4>h{4>ym4>. 


(4.11) 


Multiplying  both  sides  of  the  integral  Equation  (4.10)  by  the  testing  functions  and  using 
the  inner  product  of  Equation  (2.7)  results  in 


\%nc 

tan 


'e  wu.EWi))  • 

^  •  3  '  n 


(4.12) 


Because  of  the  Fourier  orthogonality,  the  inner  product  is  zero  for  m  ^  n  and  non-zero  for 
m  =  n.  From  now  on,  only  the  subscript  n  is  used.  Equation  (4.12)  expanded  in  matrix 
form  is 


1 

si 

'53* 

€ 
i _ 

(W2  ,EZn) 

== 

.  (Wn,E™)  . 

n 

{WuIiL(Ji))  <1 WUI2L(J2 )> 
(W2,hL(J i))  (W2,I2L(J2)) 


(WUINL(JN)) 
( W2JNL(JN )) 


[  (WnJxL^)) 

(wn,i2l(J2))  •••  (wn,inl(Jn))  Jn 

'  (WuL(Ji)) 

(Wi,L(J2))  •••  (WUL(JN))  ‘ 

'  h  ' 

<1 w2,L(J 0) 

(W2,L(J2))  •••  (W2,L(Jn)) 

h 

.  (WN,  L{J i)) 

(1 Wn,L(J2 ))  •••  (Wjv,L(Jjv)>  . 

n 

.  In  . 

[V]n  =  [Z]n[I]n. 


(4.13) 
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Equation  (4.13)  represents,  for  each  Fourier  mode,  N  equations  and  N  unknowns,  In. 
Recall  that  W{  and  Jj  both  have  t  and  4>  terms.  This  further  expands  Equation  (4.13)  into 
a  pair  of  N  equations  for  each  component.  Instead  of  doing  this,  however,  the  mathematics 
of  linear  algebra  and  the  matrix  notation,  allow  simultaneous  solutions  for  the  t  and  (j> 
component  equations  by  vertically  stacking  the  equations.  The  impedance  matrix,  [Z], 
contains  the  (W*  +  Wf  j  •  (jj  +  j/)  cross-terms, 


- 1 

1 - 1 

3$. 

_ 1 

"  [Z*] 

to 

! 

g 

_ 1 

1 - 

s-s- 

l _ 

t  N  ( 

\z**\ 

[tf] 

where, 

V*n  = 

v+  =  (Wf,Etn) 

Z%  =  {W- ,  L(Jj)) 

K*  =  (Wf,L(j})) 
z£  =  (Wf,L(Jj)) 
z£ *  =  (Wf,L(Jf)) 

ii  =  i] 

It  =  if.  (4.15) 

4  A  The  Impedance  Matrix 

Performing  the  inner  products  for  [Zn]  in  Equation  (4.15)  results  in  the  expressions 
given  in  Equations  (4.16)  through  (4.19)  below.  The  integration  in  has  been  performed 
analytically.  The  equations  are  found  to  match  those  given  in  by  Mautz  [9]  if  the  operator, 
T,  and  the  incident  field,  Eznc ,  are  scaled  by  I/17,  where  77  is  the  intrinsic  inpedance.  The 
integrations  in  ft  can  be  written  as  integrations  over  half  a  period  due  to  the  even  and  odd 
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symmetry  of  their  integrands, 


It  L  J  W  t^2  sin(v)  sinK)  +  Gi  cos(u) 


cos(u')] 


1 1  -pm  {*V/j(0<3s«ta(>0  +  ^gp  [//,(*')]  <?.}  *'* 
U/fiV)  {*V/i(‘)G»»in(»)  + MWIft}*1* 

f  f  jpfi(t)p,fj(t')  \  ^G2 - 7G1  ]  dt'dt, 

JtJtf  V  pp  J 


(4.16) 

(4.17) 

(4.18) 

(4.19) 


with  the  following  definitions, 

/*7T  p-jkR 

Gi  =  Jo  — —  cos{n4>')d<j)'  (4.20) 

pTT  p—jkR 

G2  =  /  -rD  cos {<j>')  co&{n^')d(j)'  (4.21) 

Jo  fcTt 

/“7T  p-jkR 

G3  =  -  sin(</>r)  sm{n<)>')d<j>'  (4.22) 

Jo  kR 


£ 


(/>  —  />')2  +  (z  —  27)2  +  App'sin 2 


(4.23) 


The  quantities  p,  z,  and  v  are  functions  of  arclength,  t  as  shown  in  Figure  3.1.  Similarly,  p', 
z' ,  and  1/  are  functions  of  t'.  The  functionality  in  the  expression  for  R  is  written  as  <f>' 
only  because  of  the  rotational  symmetry.  This  is  shown  by  a  simple  change  of  variables  in 
the  integrals  involving  R.  Note  that  replacing  ij  by  ji  in  Equations  (4.16)  through  (4.19) 
implies  replacing  t  dependent  quantities  by  t'  dependent  quantities,  and  vice-versa.  The 
following  symmetries  are  apparent  (a  result  of  Galerkin  testing)  and  will  reduce  the  number 
of  computations  required. 


(4.24) 
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4-2  The  Voltage  Matrix 

The  inner  products  of  Equation  (4.15)  for  the  voltage  matrix,  [Vn],  give  Equa¬ 
tions  (4.25)  and  (4.26).  The  voltage  matrix  represents  the  incident  excitation  of  the  system, 

[V*].  =  ft  kpfi(t)  (t  •  e-k-f~n*d<l>dt  (4.25) 

[*],  =  JJ  kpfi{t)  (bhinc)  e-lf-nU<t>dt,  (4.26) 

where  u%nc  is  the  polarization  of  the  incident  electric  field,  6tnc  or  0inc . 


4 >3  Scattered  Field  Equation 

The  current  density  induced  on  the  surface  of  the  body  re-radiates  to  create  the 
scattered  field.  Using  the  current  density  determined  by  the  Moment  Method  solution,  an 
equation  can  be  derived  based  on  reciprocity  which  describes  the  scattered  field.  A  row 
matrix  can  be  itentified  in  the  equation  that  parallels  the  voltage  matrix  of  the  previous 
section  allowing  it  to  be  filled  using  the  same  computer  subroutines.  For  this  it  is  useful 
to  rewrite  the  scalar  potential  in  Equation  (3.4)  in  terms  of  the  magnetic  vector  potential, 


$e=  -- T^— V-i4. 
June 


The  expression  for  Es  becomes 


(4.27) 


Es  =  -j— V(V  •  A)  -  juA.  (4.28) 

ufie 

Balanis  [1]  notes  that  the  first  term  in  Equation  (4.28)  contains  only  variations  of  the  order 
1/r2,  1/r3,  1/r4,  etc.  The  variation  of  order  1/r  is  contained  in  the  second  term  and  is  the 
dominant  variation.  In  the  far-field  (large  r),  the  first  term  is  neglected  and  the  remaining 
radial  component  becomes  negligible  compared  to  the  6  and  <f>  components. 

The  field  radiated  in  the  far-zone  is  due  to  the  current  density  on  the  surface.  By 
reciprocity,  this  is  equivalent  to  saying  that  the  field  on  the  surface  is  due  to  radiation  from 
an  infinitesimal  current  element  in  the  far-zone.  The  field  at  the  surface,  S,  is  determined 
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by  Equation  (4.28)  and  is  shown  by  Harrington  [5]  to  be 


Ers(I0u)  =  I0ue-^, 


47r  rr 


(4.29) 


where  Iqu  is  the  current  element  vector  pointing  in  the  9  or  <f>  direction,  rr  is  the  distance 
from  the  far  field  point  to  the  origin,  r  is  the  distance  from  the  origin  to  a  point  on 
the  surface,  and  kr  is  the  propagation  vector  of  the  plane  wave  coming  from  the  current 
element.  Taking  the  inner  product  of  Ers  and  the  current  density  expansion  gives  the 
magnitude  of  the  scattered  field  due  to  J  in  either  the  9  or  <f>  direction, 


Es  _  - jkrje jfc?v  y, 
47T7Y  ^ 


K 

T 

i - 

i _ 

Rt  _ 

i 

Jn<t>r 


(4.30) 


where  the  superscript  T  signifies  matrix  transpose  and  the  elements  of  the  row  matrix, 

[■Rn]  =  [Rt  Rt ]  are  siven  by 


kpfi(t)  (t  •  e  k'r+n^d<j>dt  (4-31) 

[Rt].  =  (4.32) 

where  iir  is  the  receiver  polarization  unit  vector,  9r  or  <jf .  Comparison  of  Equations  (4.31) 
and  (4.32)  with  Equations  (4.25)  and  (4.26)  show  that  when  the  measured  field  unit  vector 
is  replaced  by  the  incident  field  unit  vector,  the  only  difference  is  the  sign  on  n.  The 
following  unit  vector  dot  products  are  derived  from  the  geometry  and  are  tabulated  here 
for  use  in  the  equations  for  [Vn]  and  [i?n], 


t  •  6r't  =  —  sin(0r,<)  cos(u)  +  cos(0r’*)  sin(v)  cos(^>) 
0-  6r yt  =  —  cos(0r’*)  sin(</>) 
i  •  cjf'1  —  sin(u)  sin(<^>) 

.  iff'1  —  cos(<fi) 

—kVit  •  r  =  kz cos(0r’*)  +  kpsin(0r,t)  cos(^>), 


(4.33) 
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where  6r't  is  the  transmitted  and  received  electric  field  elevation  angle,  respectively  (Fig¬ 
ure  3.2).  The  angles  </>  and  v  are  as  shown  in  Figure  3.1. 


4-4  Three-dimensional  Scattering  Cross-section 

A  useful  measure  of  the  scattering  properties  of  a  body  is  the  scattering  cross-section . 
This  quantity  describes  the  area  that  would  be  required  to  intercept  the  amount  of  power 
that  a  source  would  have  to  scatter  isotropically  to  create  the  same  power  density  at 
the  receiver  as  does  the  actual  scattering  body  [2].  Once  the  scattered  electric  field,  E8 , 
is  determined,  this  scattering  cross-section  can  be  calculated  and  is  useful  to  compare 
the  scattering  of  bodies  of  various  sizes  and  shapes.  For  a  three-dimensional  body,  the 
scattering  cross-section,  <7,  is  defined  for  plane- wave  incidence  as  a  ratio  of  the  power 
in  the  scattered  field  to  the  power  in  the  incident  field  [2].  Noting  that  the  power  is 
proportional  to  the  field  magnitude  squared, 


cr  =  lim  47r r2 
r— yoo 


fiinc 


(4.34) 


Writing  the  expression  for  scattered  field  in  Equation  (4.30)  as  a  one-sided  sum  in 
n  using  the  even  and  odd  symmetries  of  the  real  and  imaginary  components  of  and 
allowing  <j>r  to  equal  zero  (monostatic  measurement)  gives, 


(4.35) 


Also,  as  an  aid  in  the  comparison  of  objects  of  different  sizes  at  different  frequencies,  a  can 
be  normalized  with  respect  to  the  wavelength,  A  =  2n/k,  of  the  incident  field.  Substituting 
in  \Emc\  =  hr)  from  Equation  (3.3),  the  normalized  far  field  (r  — >  oo)  cross-section  reduces 


to 


o  1  27 rrrEs(J) 

A2  47T3  7] 


(4.36) 
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4-5  Summary 

The  boundary  condition  on  the  tangential  electric  fields  at  the  surface  of  the  PEC 
body  was  used  in  Chapter  III  to  develop  an  integral  equation  to  which  the  Moment  Method 
was  applied  in  this  chapter.  Matrix  elements  were  defined  in  terms  of  an  inner  product  with 
Galerkin  testing  functions  for  a  voltage  matrix  and  an  impedance  matrix.  The  elements  of 
a  scattering  matrix  were  also  defined.  These  elements  were  computed  using  a  computer 
and  the  system  of  equations  was  solved.  The  results  are  reported  in  Chapter  V. 
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V.  Results 


Both  perfectly  electrically  conducting  spheres  and  closed  right-circular  cylinders  are  used 
as  calibration  objects  on  radar  ranges.  An  exact  solution  exists  for  the  scattering  from 
spheres.  Desire  for  the  ability  to  accurately  calculate  the  scattering  from  cylinders  is,  in 
part,  the  motivation  for  this  research.  The  equations  of  the  previous  chapter  are  derived 
for  arbitrary  body  of  revolution  geometries  and  choice  of  expansion  function.  However, 
since  the  exact  solution  is  known  for  a  sphere,  it  is  chosen  for  the  proof  of  concept  data 
presented  in  this  chapter. 

The  impedance  matrix  and  the  voltage  matrix  were  created  for  spheres  of  various 
radii  by  implementing  the  equations  of  the  previous  chapter  in  a  computer  code  written  for 
Matlab  5.  The  computer  code  was  executed  on  a  Sun  Microsystems  Sparc  20  workstation 
running  the  Solaris  2.6  operating  system.  A  function  is  called  by  the  main  program  to 
generate  the  body  of  revolution  and  compute  an  array  of  t  and  ^  current  density  expansion 
coefficients  given  input  parameters  of  sphere  radius,  maximum  sub/entire-domain  mode, 
fourier  mode  range,  backscatter  angle,  and  specification  of  sub-  or  entire-domain  method. 
The  far-field  scattering  cross-section  is  also  computed.  The  source-code  is  included  in 
Appendix  A  of  this  document  and  can  be  easily  modified  for  other  geometries  and  basis 
functions. 

In  this  chapter,  computed  results  are  presented  for  the  scattering  from  spheres  rang¬ 
ing  in  electrical  size  from  ka  =  1  to  ka  =  71.  The  scattering  is  calculated  for  subdo¬ 
main  current  expansions  using  pulse  functions  and  entire-domain  current  expansions  using 
Chebyshev  polynomials  of  the  first  kind.  Galerkin  testing  is  used  in  both  cases.  The  data 
shows  that  the  entire-domain  expansions  converge  more  accurately  and  in  fewer  modes 
than  the  subdomain  expansions.  A  comparison  showing  the  potential  computation  time 
savings  of  the  entire-domain  method  over  the  subdomain  method  is  also  presented.  It  was 
found  that  phase  approximations  used  in  the  discretization  and  numerical  solution  limit  the 
accuracy  of  the  computed  scattering  for  spheres  of  small  electrical  size.  Recommendations 
for  improving  the  computations  are  given. 

~  2tt/A  and  a  =  sphere  radius. 
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5.1  The  Computer  Code 


The  main  routine  of  the  computer  code,  MAINSPHERE,  is  written  in  modular  fash¬ 
ion  for  clarity.  The  subroutines  MAKET,  MAKEZV,  COMPRCS,  MAKELIMITS,  and 
PLOTRCS  are  written  separately  only  for  readability.  MAINSPHERE  initiates  all  vari¬ 
ables  and  calls  the  subroutines.  MAKET  creates  a  25  x  25  matrix  of  the  polynomial 
coefficients  of  the  Chebyshev  polynomial  of  the  first  kind  to  be  used  as  an  entire-domain 
basis.  MAKEZV  computes  the  impedance  matrix,  the  voltage  matrix,  and  the  measure¬ 
ment  matrix.  The  current  coefficients  are  determined  by  computing  the  inverse  of  the 
impedance  matrix  and  performing  a  matrix  multiplication.  COMPRCS  computes  the  far- 
field  scattering  cross-section  which  is  plotted  by  PLOTRCS.  MAKELIMITS  creates  an 
array  of  integration  limits  used  if  the  subdomain  method  is  chosen.  A  flowchart  of  the 
code  is  shown  in  Figures  5.1  through  5.3. 


The  integrands  of  Equations  (4.16)  through  (4.19)  for  [Z**]^.,  [z^j..,  and 

[z**]  of  the  impedance  matrix  are  written  as  functions  of  their  variables,  t,  t',  and 
<j>'.  These  functions  are  ZTTARG2,  ZTPARG2,  ZPTARG2,  and  ZPPARG2,  respectively. 
Integration  is  performed  by  the  functions  similarly  named  (without  the  suffix  ARG,  i.e., 
ZTT2,  etc.)  which  call  GQUAD3D  to  create  a  three-dimensional  array  of  integration 
points  according  to  a  Gaussian  quadrature  rule  and  evaluate  the  integrands  at  those  points. 
Matlab  allows  a  “simultaneous”  function  evalution  (2-3  sec  total  CPU  time  for  a  Z  integral 
evaluation)  of  all  points  in  the  array  using  built-in  matrix  operations.  An  early  draft  of 
the  code  attempted  to  evaluate  the  integrand  by  stepping  through  each  point  separately. 
This  took  the  same  CPU  time  per  integration  point  and  was  prohibitive.  Mautz  [8]  used 
20  integration  points  for  the  <j>'  integration  and  a  four-delta  amplitude  approximation  per 
subdomain  segment  with  satisfactory  results.  Since  MAINSPHERE  is  used  to  compute 
the  solution  using  either  subdomain  or  entire-domain  basis  functions,  Gaussian  quadrature 
is  used  to  integrate  with  respect  to  t  and  t'  as  well  as  (p .  Using  20  integration  points 
for  all  three  variables  produced  integral  convergence  to  at  least  eight  significant  digits. 
Convergence  to  more  than  ten  significant  digits  was  obtained  using  30  integration  points. 
Since  computation  time  was  comparable,  30  integration  points  for  each  variable  is  used 
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Figure  5.1  Flowchart  of  MAINSPHERE 
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MAKEZV 


Figure  5.2  Flowchart  of  MAKEZV 
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for  all  results  presented  here.  The  symmetries  expressed  in  Equation  (4.24)  are  exploited 
thereby  cutting  computation  time  roughly  in  half. 

The  integrands  of  equations  for  [V/]  and  [vf]  are  similarly  written  as  functions  of 
their  variables,  t  and  <f>.  Integration  is  performed  by  a  two-dimensional,  30  point  Gaussian 
quadrature  method.  The  measurement  vectors,  [r\\  and  [i?f]  are  produced  by  integrating 
|V/]  and  [vf\  with  the  Fourier  mode  index  negated,  —  n. 

The  function  GEO  returns  the  parameters  rho ,  zee,  and  vee  which  define  the  body  of 
revolution  generating  arc  as  a  function  of  arclength,  t  or  t '  as  appropriate.  The  parameters 
rho  and  zee  represent  the  lengths  of  the  p(t)  and  z(t)  components  of  a  vector  from  the 
origin  to  a  point  on  the  generating  arc  defined  by  t  as  shown  in  Figure  3.1.  The  parameter 
vee  represents  the  angle  between  the  rotational  axis,  £,  and  the  tangent  vector  to  the 
generating  arc  at  t ,  pointing  in  the  direction  of  increasing  t.  The  sign  of  vee  is  defined  and 
used  as  a  positive  number  if  the  tangent  at  t  points  away  from  z  and  a  negative  number  if 
the  tangent  at  t  points  toward  £.  The  parameters  that  define  a  sphere  in  positive  z  are, 


p  =  a  sin  (t/a) 


z  =  —  acos(t/a)  +  a 

.  _i  f  A p  \ 

v  =  sin  — , —  =  ■—■-==  , 

\  y/Ap2  +  Az2  ) 


(5.1) 


where  a  is  the  radius  and  A p  and  A z  are  increments  in  p  and  z  due  to  an  infitesimally 
small  increase  in  the  arclength,  t ,  at  the  point  on  the  generating  arc  specified  by  a  given  t . 

The  function  BASIS  takes  a  mode  and  an  arclength  as  inputs  and  returns  the  value 
of  the  requested  basis  function  mode  at  that  arclength.  This  is  done  via  a  mapping  of 
the  arclength  domain  in  tf  to  the  Chebyshev  basis  function  domain  of—  I<s<+1.  For 
impedance  integral  calculation,  the  functions  DERIVTP  and  DBASIS  are  used  to  com¬ 
pute  the  value  of  the  derivative  of  the  quantity  p'fj(tf).  The  companion  functions  TEST, 
DERIVT,  and  DTEST  compute  similar  quantities  with  respect  to  the  testing  functions. 
Since  Galerkin  testing  is  used,  these  routines  are  exactly  the  same  as  the  basis  function 
routines  and  can  be  combined.  The  code  was  written  to  call  them  separately,  however,  for 
readability  and  demonstration  purposes. 
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5.2  The  Currents 


The  subroutine  MAKEZV  fills  the  [Z]  matrix  and  the  [F]  matrix.  The  current  density 
expansion  coefficients  Ja  from  Equation  (4.15)  are  determined  by  the  inverse  impedance 
matrix,  [/]  =  [Z]~x  [V],  where  [I]  is  the  transpose  of  [/*  J^j .  The  elements  of  both  [I1] 
and  [7^]  are  the  coefficients  of  the  corresponding  terms  in  the  current  expansions  in  the  t 
and  <j>  directions  for  a  given  Fourier  mode,  n.  The  total  currents  are  determined  by  (using 
function  and  variable  names  from  the  computer  code),  for  each  n, 

edmax 

Jn  =  i'  Y  It{j)*basis(j,t') 
j~  0 

6(1 TTLCL  3? 

Jn  =  ft  Y  (5.2) 

i=o 

where  tf  is  a  vector  of  samples  across  the  arclength,  0  <  tf  <  tmax.  Edmax  is  the  maximum 
number  of  entire-domain  modes  used.  Adding  across  all  n  gives, 

fmax 

J1' =  t'  Y  *exp(j  *n*<j>'),  (5.3) 

n— 0 

where  fmax  is  the  maximum  number  of  Fourier  modes.  It  is  convenient  to  evaluate  the 
currents  along  the  <£'  =  0  longitude  (in  the  XZ- plane),  because  the  exponential  becomes 
unity.  Figure  5.4  shows  the  <£  directed  current  densities  represented  by  each  entire-domain 
mode  for  a  sphere  of  ka  —  4.0  excited  by  a  ^  directed  incident  electric  field  when  n  —  2 
evaluated  at  broadside  (9  =  k/2).  There  is  a  similar  set  of  currents  for  every  Fourier 
mode  used  in  the  current  expansion.  The  top  two  subplots  show  the  contribution  to  the 
current  density  of  each  of  the  T0-T5  entire-domain  modes.  The  modes  are  recognizable 
as  the  function  of  Figure  3.3  scaled  by  the  corresponding  element  of  [/*].  The  scaling  for 
some  modes  is  very  small  or  zero,  minimizing  their  effect  on  the  overall  current.  Only  the 
modes  that  serve  to  accurately  represent  the  functionality  of  the  actual  current  density 
as  described  by  the  integral  equation  are  selected  by  the  method.  The  bottom  subplot  of 
Figure  5.4  shows  the  sum  of  each  entire-domain  mode  above.  The  solid  line  here  represents 
the  current  density  magnitude  for  the  n  =  2  Fourier  mode.  The  total  current  density  is  the 
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^-Component  of  Current  Density  for  0=90  deg 
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Total  current  density  over  all  ed/f  modes,  9=90  deg,  ka=4 


Figure  5.5  Total  i  and  <j)  directed  current  densities  (n  =  0-5)  for  ka  =  4.0  illuminated 
by  Efc  at  6  =  90° 

sum  of  all  entire-domain  contributions  for  all  Fourier  modes.  The  total  <j>  directed  current 
density  for  a  sphere  of  ka  =  4.0  evaluated  at  broadside  (0  =  7r/2)  is  shown  in  Figure  5.5. 
A  sphere  of  radius  a  =  4.0 /k  has  a  circumference  of  four  wavelengths.  The  current  density 
around  the  sphere  exhibits  more  oscillations  than  does  a  sphere  of  smaller  radius.  Consider 
the  smoother  total  current  density  for  a  sphere  of  radius  a  =  3.1/A:  as  shown  in  Figure  5.6. 
The  circumference  is  now  approximately  three  wavelengths.  More  Fourier  series  expansion 
terms  of  higher  oscillation  orders,  are  required  to  accurately  describe  the  current 

density  for  spheres  of  larger  electrical  size. 

The  (j>  directed  current  density  coefficients  for  each  entire-domain  mode  are  plotted 
in  Figure  5.7  for  sphere  sizes  of  ka  =  0.9,  ka  =  2.1,  and  ka  =  4.0.  The  coefficients  are 
grouped  by  Fourier  mode  in  sequential  order  of  Inj  =  Jno  through  Iu\q.  The  vertical  axis  of 
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Total  current  density  over  all  ed/f  modes, 9=90  deg,  ka=3.1 


Figure  5.6  Total  (j)  directed  current  density  (n  =  0-5)  for  ka  =  3.1  illuminated  by  E™c 
at  6  =  90° 
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4>—Current  Density  Coefficient 


Fourier  mode,  n 

Figure  5.7  <f)  directed  current  density  coefficients  for  the  Ib-Tio  entire-domain  modes, 

grouped  by  Fourier  mode 

each  plot  is  scaled  to  show  the  appropriate  detail  for  each  sphere  size.  It  is  easily  seen  that 
as  the  sphere  size  increases,  more  Fourier  modes  are  required  to  capture  the  total  energy 
of  the  current  density.  In  other  words,  more  Fourier  modes  must  be  considered  before  the 
current  density  coefficients  decrease  sufficiently  in  amplitude.  Also,  for  a  given  Fourier 
mode,  it  is  obvious  that  most  of  the  energy  is  captured  in  the  To  through  T5  entire-domain 
modes. 

Figure  5.8  shows  the  t  directed  current  densities  represented  by  each  entire-domain 
mode  for  a  sphere  of  ka  =  4.0,  n  =  2  when  excited  at  broadside  ( 6  =  ir/2)  by  a  <j>  polarized 
incident  electric  field.  There  is  a  similar  set  of  currents  for  every  Fourier  mode  used  in 
the  current  expansion.  Note  that  as  in  the  directed  current,  some  entire-domain  modes 
are  emphasized  while  others  are  minimized.  As  before,  this  serves  to  represent  the  actual 
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t-Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=2,  E^nc 
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t— Current  Density  Coefficient 


Fourier  mode,  n 


Figure  5.9  t  directed  current  density  coefficients  for  the  To-Tio  entire-domain  modes, 
grouped  by  Fourier  mode 

current  density  functionality  in  this  direction.  The  sum  of  the  entire-domain  modes  is 
represented  in  the  bottom  plot  of  Figure  5.8  where  the  functionality  is  seen  to  be  near 
zero  across  most  of  the  arclength.  The  total  t  directed  current  density  for  the  sphere  of 
ka  =  4.0  evaluated  at  broadside  (0  =  tt/2)  is  shown  in  Figure  5.5.  The  t  directed  current 
density  coefficients  for  each  entire-domain  mode  are  plotted  in  Figure  5.9  for  sphere  sizes 
of  ka  =  0.9,  ka  =  2.1,  and  ka  =  4.0.  The  trends  are  similar  to  those  of  the  cj)  directed 
current  density  coefficients.  The  current  densities  for  all  modes  of  the  ka  =  4.0  sphere  in 
both  the  t  and  <j>  directions  are  given  in  Appendix  B. 

Wood’s  subdomain  analysis  of  the  current  density  excited  on  the  surface  of  a  spherical 

A  A 

cavity  finds  that,  with  the  exception  of  the  n  =  1  Fourier  mode,  both  the  t  and  6  current 
density  vanishes  at  the  poles  [16].  The  amplitude  coefficient  for  each  term  in  the  current 
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density  expansion  is  determined  by  the  Moment  Method  computation  based  on  the  need 
for  the  frequency  content  that  term  offers.  Emphasis  of  a  particular  term  is  based  on  a 
weighted  average  over  the  domain  of  the  portion  of  the  generating  arc  to  which  the  term 
applies.  In  a  subdomain  method,  the  coefficient  for  the  function  which  represents  the 
current  density  amplitude  at  the  poles  can  be  zero  without  respect  to  any  other  segment 
of  the  generating  arc.  In  an  entire-domain  method,  the  domain  of  each  mode  is  the 
entire  length  of  the  generating  arc,  and  the  amplitude  coefficient  of  each  mode  is  the  best 
choice  across  the  domain  of  the  arc.  In  choosing  Chebyshev  expansion  functions,  there  is  a 
potential  disadvantage  in  that  each  mode  has  a  positive  or  negative  finite  amplitude  at  both 
ends.  The  only  way  for  the  current  to  vanish  at  the  poles  is  by  cumulative  addition  of  many 
modes  scaled  by  their  coefficients.  This  is  not  a  disadvantage  of  using  entire-domain  basis 
functions,  but  rather,  an  indication  that  another  entire-domain  family  of  functions  may  be 
more  appropriate  for  this  specific  geometry.  Given  enough  Chebyshev  modes,  however,  it 
is  expected  that  the  current  would  approach  the  anticipated  value  across  the  entire  domain 
of  the  sphere,  including  at  the  poles.  Because  of  the  decreased  amplitude  of  the  higher 
order  modes,  many  more  may  be  required  than  that  necessary  for  determining  the  current 
density  at  points  away  from  the  poles.  Had  the  fact  that  the  current  should  vanish  at  the 
poles  been  considered  at  the  beginning  of  this  research,  a  Fourier  expansion  (exponential 
basis  functions)  may  have  been  chosen  for  the  expansion  of  the  current  density  along  the 
arclength  as  well  as  in  <j>.  This  would  have  allowed  a  build-up  of  the  needed  frequency 
content  satisfied  by  the  lower  order  modes  without  adversely  affecting  the  current  density 
at  the  poles. 

When  the  sphere  is  excited  by  a  plane  wave  incident  at  an  angle  other  than  broadside, 
the  poles  shift  and  are  no  longer  at  the  ends  of  the  generating  arc.  A  corresponding  shift 
in  the  current  density  distribution  along  the  arclength  is  expected.  Figure  5.10  shows 
the  accumulation  of  current  density  for  ka  =  4.0  and  E™c  incident  at  0  =  7t/3.  Each 
entire-domain  mode  for  the  n  —  2  Fourier  mode  of  the  <f>  component  of  current  density 
is  plotted  in  the  top  two  subplots.  A  different  weighting  of  the  entire-domain  modes 
is  seen  when  compared  to  the  broadside  case  of  Figure  5.4.  This  weighting  serves  to 
enhance  the  current  density  on  the  generating  arc  near  the  ^  —  7t/3  region  and  contains 
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Total  current  density  over  all  ed/f  modes, 0=60  deg,  ka=4 


Figure  5.11  Total  t  and  (j)  directed  current  densities  (n  =  0-5)  for  ka  =  4.0  illuminated 
by  Efc  at  9  =  60° 

the  appropriate  frequency  content  for  this  excitation.  Note  that  there  is  still  a  rise  in  the 
current  density  near  the  ends  of  the  arclength  due  to  the  finite  sum  of  the  Chebyshev 
functions.  The  cumulative  total  over  all  entire-domain  and  Fourier  modes  (5  modes  in 
each  case)  is  plotted  in  Figure  5.11.  The  expected  shift  in  the  pattern  is  seen  with  a  peak 
at  t  —  1.33  corresponding  to  the  illumination  angle,  9 . 

5.3  The  Scattering  Cross-Section 

5.3.1  Assessment  of  Accuracy.  The  solution  to  the  integral  equations  given  in 
Chapter  IV  is  the  vector  [I]  which  contains  the  coefficients  to  the  expansion  of  the  cur¬ 
rent  density,  J.  The  accuracy  of  the  solution,  therefore,  is  best  measured  by  comparison  to 
known  currents.  Production  codes  (subdomain  basis)  such  as  CICERO  and  JRMBOR  that 
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are  currently  in  use  calculate  the  current  internally  and  output  only  the  electromagnetic 
field  components  or  the  scattering  cross-section  as  the  quantity  of  interest.  Although  the 
currents  computed  in  these  codes  are  also  an  approximation  limited  by  the  smallness  of  the 
subdomain  discretization  used,  they  would  serve  as  a  useful  check  of  what  is  obtained  here. 
Without  modification  to  the  codes  themselves,  however,  the  currents  are  not  available  to 
the  user.  Mautz  [9]  assesses  accuracy  by  defining  a  metric  that  relates  the  currents  result¬ 
ing  from  his  magnetic-  and  electric-field  integral  equations  to  his  combined-field  integral 
equation  (a  linear  combination  of  the  two).  His  purpose  is  to  show  a  deviation  from  an 
assumed  norm  of  stray  data  points  resulting  from  physical  resonances  of  the  body.  This 
is  suitable  for  showing  how  well  the  combined-field  integral  equation  discriminates  against 
the  resonances. 

Accuracy  in  the  computation  of  the  current  is  paramount  to  obtaining  an  accurate 
scattering  cross-section.  The  averaging  effect  of  the  scattering  cross-section  calculation, 
however,  allows  for  some  variation  in  current  quantities  without  largely  affecting  the  scat¬ 
tering  cross-section  itself.  It  will  here  suffice  to  use  the  far-field  scattering  cross-section  as 
the  figure  of  merit  and  compare  that  computed  by  MAINSPHERE  using  both  subdomain 
and  entire-domain  expansions  of  the  current  density  to  that  produced  by  the  Mie  series  [5] 
which  is  an  exact  solution. 

5.3,2  Subdomain  Data.  The  code  MAINSPHERE  and  its  supporting  subroutines 
and  functions  listed  in  Appendix  A  were  used  to  compute  the  scattering  cross-section  of  a 
perfectly  electrically  conducting  sphere  of  electrical  size  ranging  from  ka  =  0.5  to  ka  =  5.0. 
Approximately  30  data  points  were  computed  in  this  range.  For  this  formulation,  the 
surface  current  density  was  expanded  in  the  tf  and  </>'  directions  using  unit  amplitude 
pulses  as  basis  functions.  The  current  was  expanded  additionally  by  a  common  Fourier 
series  exponential  in  <j).  MAINSPHERE  was  run  in  the  subdomain  (‘ sd ’)  mode.  The 
resulting  quantities  are  normalized  with  respect  to  electrical  size  and  wavelength  so  as  to 
oscillate  about  a  constant  value.  Comparison  is  made  to  the  exact  normalized  scattering 
cross-section  calculated  by  means  of  the  Mie  series  given  by  Harrington  [5]. 
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Backscatter  from  PEC  Sphere 


Figure  5.12  Backscatter  from  a  PEC  sphere  using  5  subdomain  segments,  illuminated 
at  9  =  90°  by  Efc. 

The  first  data  set  shows  the  scattering  cross-section  computed  with  five  subdomain 
segments  on  the  generating  arc  for  the  selected  range  of  electrical  sizes.  Five  Fourier  modes 
were  used  in  the  <j>  expansion  for  this  and  all  other  data  presented.  The  data  points  are 
shown  in  Figure  5.12  as  discrete  points  (circles).  The  data  points  are  connected  by  a 
dotted  line  to  make  the  sequence  more  apparent  where  neighboring  points  follow  the  same 
pattern.  The  line  is  not,  however,  intended  to  indicate  an  accurate  interpolation  between 
points.  An  occasional  stray  point  gives  the  impression  that  addition  of  another  point  near 
it  would  result  in  nearly  the  same  magnitude,  but  this  is  not  necessarily  the  case.  Although 
data  points  in  the  ka  =  2.25  to  ka  =  2.75  region  approach  the  exact  solution  and  points 
in  the  ka  =  3.5  to  ka  =  4.0  approach  the  exact  solution  to  a  lesser  degree,  the  data  overall 
is  somewhat  erratic.  It  is  clearly  obvious  that  five  subdomain  segments  are  insufficient 
to  accurately  represent  scattering  cross  section.  A  large  number  of  data  points  differ  by 
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Figure  5.13  Backscatter  from  a  PEC  sphere  using  5  entire-domain  modes,  illuminated 
at  6  =  90°  by  Efc . 

=L3-10dB.  The  largest  deviation  is  a  cross-section  of  0.001  at  ka  =  3.1  which  differs  from 
the  Mie  series  by  -27dB.  Mautz  shows  that  at  least  fifteen  subdomain  sections  are  required 
to  model  the  current  density  with  accuracy  sufficient  to  match  the  Mie  series  [9]. 

5.3.3  Entire-Domain  Data .  The  next  data  set  shows  the  scattering  cross-section 
computed  with  five  entire-domain  modes  on  the  generating  arc  for  the  selected  range 
of  electrical  sizes.  The  data  points  are  shown  as  circles  connected  by  a  dotted  line  in 
Figure  5.13.  It  is  readily  apparent  that  five  entire-domain  modes  result  in  a  much  closer 
conformity  to  the  Mie  series  than  do  an  equivalent  number  of  subdomain  segments.  The 
entire-domain  basis  functions  are  better  able  to  accurately  represent  the  smoothness  and 
functionality  of  the  surface  current  density.  The  data  points  in  the  ka  =  3.0  to  ka  =  5.5 
region  closely  follow  the  Mie  series  curve  with  less  than  0.5dB  error.  The  region  from 
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ka  =  1.75  to  ka  =  3.0  exhibits  an  error  of  2dB  to  3dB.  Electrical  sizes  smaller  than 
ka  =  1.75  differ  by  an  almost  constant  -4dB. 

At  larger  electrical  sizes  (ka  >  5.5)  the  data  begins  to  fall  away  from  the  Mie  series 
curve.  This  is  due  to  the  fact  that  as  the  sphere  gets  larger,  the  current  around  the  surface 
becomes  more  varied  and,  as  a  result,  requires  more  and  more  Fourier  expansion  modes 
(higher  values  of  n )  to  accurately  represent  the  current  oscillations  (see  Figure  5.7).  The 
current  density  at  smaller  electrical  sizes  is  more  constant  around  the  sphere  circumference, 
and  is  adequately  covered  by  a  few  Fourier  modes  (2  or  3  for  ka  <  2)  [9] . 

The  inaccuracy  in  normalized  cross-section  seen  in  this  data  for  smaller  ka  is  intro¬ 
duced  by  an  approximation  made  in  the  integrands  of  MAINSPHERE  to  allow  integration 
through  the  singularity  which  occurs  when  the  difference  between  the  radial  distance  be¬ 
tween  the  far-field  observation  point  and  the  testing/basis  functions,  R  in  the  equations, 
goes  to  zero.  In  subdomain  methods,  the  only  time  this  occurs  is  when  the  testing  func¬ 
tion  coincides  with  the  basis  function.  When  entire-domain  basis  functions  are  used,  the 
domain  of  every  mode  coincides  with  every  other.  The  value  of  R  is  zero  when  t  equals 
tf  and  <f>f  equals  zero.  It  is  also  zero  at  any  value  of  <f>f  when  t  and  tf  are  both  at  their 
minimum  or  maximum  value.  This  behavior  of  R  is  shown  in  Figure  5.14.  The  radial 
projection  in  the  Ay-plane  of  the  vector  to  a  point  on  the  generating  arc,  p,  also  goes  to 
zero  at  the  ends  of  the  generating  arc  causing  a  singular  behavior  in  the  impedance  matrix 
integrands.  The  effect  of  the  p  related  singularity  is  lessened  by  the  fact  that  the  Gaussian 
quadrature  integration  does  not  involve  a  function  evaluation  of  the  integrand  exactly  at 
the  endpoints.  Also,  the  higher  density  of  integration  points  at  the  ends  of  the  integration 
domain  allow  a  more  accurate  integration  in  those  regions.  To  keep  R  from  going  to  zero 
in  the  integrations,  an  incremental  distance,  skosh ,  is  added  to  R.  This  method  is  called 
the  method  of  equivalent  separation  and  has  been  used  in  scattering  problems  [9].  The 
distance  R  is  kept  slightly  above  zero  allowing  integration  without  extraction  and  separate 
handling  of  the  singularity.  However,  since  R  represents  a  path  of  phase  accumulation, 
a  phase  error  is  introduced  into  the  integrands.  The  smaller  the  sphere  gets  in  electrical 
size,  the  larger  this  phase  error  gets  relative  to  the  sphere  size.  It  was  found  that  using 
this  equivalent  separation  for  spheres  of  ka  <  2,  the  integrals  converge  to  an  approximate 
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Amplitude  Amplitude 


Figure  5.14  Value  of  R  as  a  function  of  t,  t',  and  <f>' . 
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INTEGRAL  OF  ZTT2  VS.  Rmin 
for:  ka=4,  n=2,  test/basis=3,  variable  #  integration  points 


Figure  5.15  The  integral  of  Zu  as  a  function  of  minimum  R  for  various  orders  of  Gaussian 
quadratrue  integration. 

solution,  but  not  to  the  Mie  series  curve.  Re-computing  with  additional  entire-domain  or 
Fourier  modes  did  not  improve  the  results  in  this  region.  Reducing  the  minimum  phase 
distance,  skosh ,  brought  the  cross-section  points  slightly  closer  to  the  Mie  series  curve,  but 
did  not  solve  the  problem.  Also,  increasing  the  order  of  Gaussian  quadrature  quickly  be¬ 
came  computation  time  inhibitive.  Figure  5.15  is  a  plot  of  the  integrated  amplitude  of  Zu 
versus  minimum  phase  distance  for  a  sphere  of  ka  =  4.0,  the  n  =  2  Fourier  mode,  and  the 
Ts  testing  and  basis  modes.  Data  is  plotted  for  20,  30,  40,  and  50  point  Gaussian  quadra¬ 
ture  integration.  As  the  Gaussian  quadrature  order  is  increased,  there  is  less  variation  in 
number  given  by  the  integration.  This  is  seen  by  the  curves  getting  closer  together  as  the 
order  is  increased.  As  the  minimum  R  is  varies,  the  integration  yields  a  result  that  does 
not  converge.  A  similar  behavior  is  also  seen  when  the  testing  function  and  basis  function 
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are  not  the  same  although  ,in  general,  the  variation  is  less  pronounced.  It  is  concluded 
that  accuracy  is  only  improved  by  a  more  rigorous  handling  of  the  singularity  by  means 
of  singularity  extraction  and  separate  evaluation  of  the  integrals  in  the  region  surrounding 
the  singularity. 

5-4  Computation  and  Matrix  Fill  Time 

It  is  seen  that  increasing  the  number  of  modes  used  to  compute  the  system  of  ma¬ 
trices  results  in  a  more  accurate  solution  for  the  unknown  current  and,  therefore,  a  more 
exact  calculation  of  scattering  cross-section.  A  major  advantage  of  using  entire-domain 
basis  functions  to  expand  the  current  density  is  that  additional  modes  can  be  added  by 
superposition  without  complete  re-computation  of  each  matrix  element.  The  integration 
limits  on  the  inner  product  for  the  entire-domain  method  remain  the  same  for  every  com¬ 
bination  of  testing  and  expansion  mode,  i  and  j.  In  other  words,  each  mode  has  as  its 
domain,  the  entire  length  of  the  generating  arc.  Higher  order  modes  are  simply  added  in 
parallel  to  lower  order  modes.  This  is  not  true  for  subdomain  basis  functions.  The  inte¬ 
gration  limits  on  the  inner  product  for  the  subdomain  method  correspond  to  the  length 
of  each  segment.  Each  segmental  subdomain  mode  adds  to  the  others  in  series .  In  order 
to  increase  the  number  of  segments  in  the  subdomain  discretization,  the  length  of  every 
segment  must  change.  This  requires  complete  re-computation  of  all  matrices.  As  a  result, 
increasing  the  order  of  the  discretization  to  improve  the  accuracy  of  the  solution  wastes 
all  computation  time  up  to  that  point.  A  comparison  of  matrix  fill  times  and  computation 
times  for  the  entire-domain  method  and  the  subdomain  method  is  shown  below. 

Elapsed  computer  processor  time  was  tracked  through  all  trials  for  the  following 
computations:  [Z|j],  [Zjf],  [. zf? *],  [V*],  \VP],  and  [Z]  inversion.  The  time  to  add  an  addi¬ 
tional  mode  was  also  tracked.  These  times  are  shown  in  Table  5.1.  The  majority  of  the 
total  computation  time  is  in  the  evaluation  of  the  integrals  to  determine  the  elements  of 
the  impedance  matrix,  \Z\.  Processor  time  averages  3.0  seconds  per  integration.  Slightly 
more  than  half  of  the  elements  of  [Z\  must  be  computed  to  fill  [Z\  because  of  symmetries. 
Because  the  fill  time  for  [Z]  dominates,  proper  choice  of  basis  function  (entire-domain  vs. 
subdomain  and  actual  function)  offers  an  important  advantage. 
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Element  ! 

Time 

a,  A-, 

3.5  sec 

3.0  sec 

2.5  sec 

0.5  sec 

0.5  sec 
0.5  sec 

Table  5.1  Average  integration  times 


Based  on  the  computation  time  observations,  the  additional  time  required  to  add 
a  mode  to  the  impedance  matrix  given  an  initial  number  of  modes  is  determined  and  is 
shown  in  Tables  5.2  and  5.3. 


Consider  an  example  of  overall  computation  time  to  fill  [Z] .  Starting  with  an  initial 
choice  of  five  entire-domain  expansion  modes,  2.75  minutes  (165  seconds)  are  required  to 
compute  the  N(2N  +  1)  =  55  symmetric  elements  of  \Z].  Suppose  that  after  determining 
the  current  and  checking  it  against  a  specified  tolerance,  more  modes  are  added.  This  is 
repeated  until  ten  modes  are  used  to  achieve  an  accurate  solution.  Every  time  a  mode  is 
added,  a  new  row  is  added  to  [Zn]  and  [Z^]  and  a  row  and  a  column  are  added  to  [Z1^] 
as  shown  in  Equation  (5.4)  for  a  total  of  4 N  +  3  new  computations,  where  N  is  the  initial 
number  of  modes.  The  x’s  represent  elements  of  previously  calculated  modes,  o’s  represent 
new  elements  that  must  be  computed,  and  #’s  represent  new  elements  filled  by  symmetry. 
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(5.4) 


This  additional  computation  time  is  added  to  the  initial  matrix  fill  time  to  get  the  overall 
computation  time.  Doubling  the  number  of  entire-domain  modes  to  ten  adds  7.75  minutes 
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ENTIRE-DOMAIN 

Initial 

Elements 

Fill 

Elements  req’d 

Add 

Total 

Computation 

modes 

to  fill 

time 

to  add  mode 

time 

CPU  time 

overhead 

1 

3 

9 

7 

21 

30 

0 

2 

10 

30 

11 

33 

63 

0 

3 

21 

63 

15 

45 

108 

0 

4 

36 

108 

19 

57 

165 

0 

5 

55 

165 

23 

69 

234 

0 

6 

78 

234 

27 

81 

315 

0 

7 

105 

315 

31 

93 

408 

0 

8 

136 

408 

35 

105 

513 

0 

9 

171 

513 

39 

117 

630 

0 

10 

210 

630 

43 

129 

759 

0 

11 

253 

759 

47 

141 

900 

0 

12 

300 

900 

51 

153 

1053 

0 

13 

351 

1053 

55 

165 

1218 

0 

14 

406 

1218 

59 

177 

1395 

0 

15 

465 

1395 

63 

189 

1584 

0 

16 

528 

1584 

67 

201 

1785 

0 

17 

595 

1785 

71 

213 

1998 

0 

18 

666 

1998 

75 

225 

2223 

0 

19 

741 

2223 

79 

237 

2460 

0 

20 

820 

2460 

83 

249 

2709 

0 

21 

903 

2709 

87 

261 

2970 

0 

22 

990 

2970 

91 

273 

3243 

0 

23 

1081 

3243 

95 

285 

3528 

0 

24 

1176 

3528 

99 

297 

3825 

0 

25 

1275 

3825 

103 

309 

4134 

0 

Based  on  average  integration  time  of  3  sec. 

All  quantities  in  seconds  unless  specified  otherwise. 

Table  5.2  Additional  time  required  per  Fourier  mode  to  add  an  entire-domain  mode 
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SUBDOMAIN 

Initial 

Elements 

Fill 

Elements  req’d 

Add 

Total 

Computation 

modes 

to  fill 

time 

to  add  mode 

time 

CPU  time 

overhead 

i 

3 

9 

10 

30 

39 

9 

sec 

2 

10 

30 

21 

63 

93 

30 

sec 

3 

21 

63 

36 

108 

171 

1.05 

min 

4 

36 

108 

55 

165 

273 

1.8 

min 

5 

55 

165 

78 

234 

399 

2.75 

min 

6 

78 

234 

105 

315 

549 

3.9 

min 

7 

105 

315 

136 

408 

723 

5.25 

min 

8 

136 

408 

171 

513 

921 

6.8 

min 

9 

171 

513 

210 

630 

1143 

8.55 

min 

10 

210 

630 

253 

759 

1389 

10.5 

min 

11 

253 

759 

300 

900 

1659 

12.65 

min 

12 

300 

900 

351 

1053 

1953 

15 

min 

13 

351 

1053 

406 

1218 

2271 

17.55 

min 

14 

406 

1218 

465 

1395 

2613 

20.3 

min 

15 

465 

1395 

528 

1584 

2979 

23.25 

min 

16 

528 

1584 

595 

1785 

3369 

26.4 

min 

17 

595 

1785 

666 

1998 

3783 

29.75 

min 

18 

666 

1998 

741 

2223 

4221 

33.3 

min 

19 

741 

2223 

820 

2460 

4683 

37.05 

min 

20 

820 

2460 

903 

2709 

5169 

41 

min 

21 

903 

2709 

990 

2970 

5679 

45.15 

min 

22 

990 

2970 

1081 

3243 

6213 

49.5 

min 

23 

1081 

3243 

1176 

3528 

6771 

54.05 

min 

24 

1176 

3528 

1275 

3825 

7353 

58.8 

min 

25 

1275 

3825 

1378 

4134 

7959 

1.0625 

hr 

Based  on  average  integration  time  of  3  sec. 

All  quantities  in  seconds  unless  specified  otherwise. 

Table  5.3  Additional  time  required  per  Fourier  mode  to  add  a  subdomain  segment 
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(465  seconds)  to  the  fill  time  for  a  total  of  10  minutes  (630  seconds), 


CPU  time  —  (I65)s  initial  modes  +  69  +  81  +  93  +  105  +  117 

=  630  seconds  to  compute  10  modes .  (5.5) 

The  same  calculation  is  made  for  an  initial  choice  of  five  subdomain  segments.  It 
takes  the  same  2.75  minutes  (165  seconds)  to  fill  the  subdomain  [Z]  as  it  does  for  entire- 
domain.  Doubling  the  number  of  segments  to  ten  adds  35  minutes  (2100  seconds)  to  the 
fill  time  for  a  total  of  37.75  minutes  (2265  seconds).  This  is  387%  more  time  to  fill, 

CPU  time  —  (I65)s  initial  segments  +  234  +  315  +  408  +  513  +  630 

=  2265  seconds  to  compute  10  segments .  (5.6) 

A  more  realistic  situation,  however,  is  that  because  of  the  smoothness  and  appropriate 
functionality  of  the  entire-domain  basis  functions,  many  more  subdomain  segments  will  be 
required  to  achieve  the  same  accuracy.  The  processing  time  for  ten  entire-domain  modes 
is  10.5  minutes  (630  seconds),  whether  ten  are  chosen  initially  or  additional  modes  are 
added  to  make  ten.  Since  ten  entire-domain  modes  achieves  the  accuracy  desired,  let  the 
subdomain  comparison  start  with  ten  segments.  Suppose  that  to  achieve  accuracy  in  this 
case,  the  number  of  segments  is  incrementally  doubled  to  twenty.  This  results  in  4.27 
additional  hours  (15,375  seconds)  for  a  total  of  4.45  hours  (16,005  seconds)  to  fill  \Z\.  This 
is  3.76  hours  more  (551%  longer)  than  if  twenty  segments  were  initially  chosen  and  4.32 
hours  more  (2440%  longer)  than  it  took  to  achieve  the  same  results  as  in  the  entire-domain 
case!  The  time  difference  is  even  worse  than  this  because  it  is  also  compounded  by  the 
number  of  Fourier  modes, 

CPU  time  =  (630)io  +  759  +  900  +  1053  +  1218  +  1395  +  1584  +  1785  +  1998  +  2223  +  2460 
=  16005  seconds  to  compute  20  segments .  (5.7) 
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A  summary  of  the  savings  in  computation  time  for  the  above  examples  is  presented 
in  Table  5.4. 


Add  time 

Total  time 

SUBDOMAIN 

Initial  modes  =  5 
Final  modes  =10 

35  min 

Initial  modes  =10 
Final  modes  =  20 

4.27  hrs 

10.50  min 
4.45  hrs 

ENTIRE-DOMAIN 

Savings  over  subdomain 

Initial  modes  =  5 
Final  modes  =10 

7.75  min 

2.75  min 
10.5  min 

4-32  hours  saved 

Table  5.4  Summary  of  CPU  time  savings  per  Fourier  mode 
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VI.  Conclusions  and  Recommendations 


This  chapter  discusses  the  achievements  of  this  research.  The  impacts  of  this  work  are 
summarized  and  recommendations  are  made  for  future  work  in  this  area. 

6.1  Conclusions 

The  primary  motivation  for  this  research  involved  improving  the  numerical  methods 
used  to  obtain  the  electromagnetic  scattering  from  perfectly  electrically  conducting  bodies 
of  revolution.  Such  bodies  are  used  to  calibrate  the  measurement  systems  of  radar  ranges 
using  a  background  subtraction  method.  Current  approaches  use  Method  of  Moments 
based  computer  codes  which  employ  subdomain  basis  expansion  functions  to  calculate 
an  approximate  solution  to  the  current  density  induced  on  the  surface  of  the  body  by 
an  incoming  plane-wave.  The  accuracy  of  the  approximation  is  limited  by  the  number 
of  subdomain  segments  used  in  the  expansion.  A  large  number  of  these  segments  are 
required  to  accurately  represent  the  current  density.  Computation  time  becomes  inhibitive 
for  an  accurate  solution,  particularly  for  a  body  of  large  electrical  size.  Incrementally 
increasing  the  number  of  subdomain  segments  used  in  order  to  increase  accuracy  requires 
re-discretization  and  loss  of  previously  calculated  results.  Entire-domain  basis  functions 
can  solve  both  of  these  problems.  Although  many  texts  make  reference  to  the  possibility 
of  using  entire-domain  basis  functions  for  expansion  of  the  current  density,  they  all  default 
to  the  apparent  simplicity  of  subdomain  basis  functions.  A  search  of  the  literature  also 
yields  few  examples  of  the  use  of  entire-domain  basis  functions. 

As  a  proof  of  the  concept  of  using  entire-domain  basis  functions  to  compute  the 
scattering  of  a  body  of  revolution,  the  unknown  current  density  on  a  perfectly  electrically 
conducting  sphere  was  expanded  using  Chebyshev  polynomials  of  the  first  kind.  The 
Method  of  Moments  was  then  used  to  solve  for  the  current  density.  Galerkin’s  method 
of  testing  allowed  an  approximate  50%  reduction  in  required  computations.  A  sphere 
was  chosen  to  enable  comparison  to  the  Mie  series,  a  known  exact  solution.  The  entire- 
domain  solution  using  the  Chebyshev  basis  was  shown  to  be  a  better  approximation  of  the 
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scattering  cross-section  in  many  fewer  modes  (about  5)  than  an  equivalent  trial  using  a 
subdomain  pulse  basis. 

Since  the  current  density  is  a  smooth  variation  over  the  surface,  the  Chebyshev 
basis  functions  were  a  better  approximation  than  discontinuous  subdomain  pulses.  This 
was  a  significant  factor  in  reducing  the  number  of  modes  required  in  the  approximation 
and  the  matrix  size  of  the  solution.  Computation  time  was  dramatically  decreased.  The 
computational  time  required  to  obtain  an  initial  solution  and  then  double  the  number  of 
entire-domain  modes  used  in  order  to  increase  the  accuracy  of  the  solution  was  found  to 
be  the  same  as  if  the  final  number  of  modes  were  used  initially.  Using  subdomain  basis 
segments  took  hours  longer — clearly  an  advantage  of  entire-domain  basis  functions  that 
overcomes  the  initial  complexity  of  their  implementation. 

6.2  Recommendations  for  Future  Research 

Although  the  results  presented  in  this  report  accomplish  the  goal  of  the  research 
effort  by  showing  the  obvious  computational  advantages  and  time  savings  of  using  entire- 
domain  basis  functions  for  the  given  problem,  there  are  several  areas  in  which  further  work 
could  follow. 

First,  an  accurate  solution  was  not  obtained  for  spheres  of  smaller  electrical  sizes. 
This  was  due  to  the  handling  of  the  weak  singular  nature  of  the  three-dimensional  free-space 
Green’s  function  by  means  of  imposing  an  equivalent  separation  in  the  minimum  distance 
between  the  testing  functions  and  the  current  expansion  functions — an  approach  that  has 
been  successful  in  other  similar  problems.  The  small  phase  approximation  introduced  was 
found  to  be  unacceptable  for  spheres  of  smaller  electrical  sizes.  Methods  exist  for  extracting 
the  singularity  and  evaluating  the  integrals  in  the  region  surrounding  the  singular  point 
separately.  Future  research  would  use  such  a  method  to  handle  the  singular  behavior  of 
the  integrands  more  rigorously.  The  Chebyshev  weighting  factor  for  true  orthogonality 
could  also  be  implemented. 

Second,  a  sphere  was  chosen  in  this  effort  to  demonstrate  the  advantages  of  entire- 
domain  basis  functions  because  an  exact  solution  exists  for  comparison.  Real  interest 
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in  this  method  lies  in  the  ability  to  accurately  calculate  the  scattering  from  bodies  of 
revolution  to  which  no  exact  solution  is  known.  Applying  appropriately  chosen  entire- 
domain  basis  functions  to  the  solution  of  scattering  from  more  complex  and  advantageous 
geometries  would  be  useful.  Examples  of  such  geometries  might  be  elliptical  bodies  of 
low  cross-section,  the  right-circular  cylinder  used  in  cylindrical  calibration  methods  for 
radar  ranges,  and  other  bodies  of  revolution  with  discontinuities  (such  as  corners)  in  their 
generating  arc.  Investigation  might  be  made  into  the  appropriateness  of  the  method  for 
bodies  with  concave  and/or  convex  regions.  Additionally,  it  would  be  useful  to  study  the 
effects  of  using  various  entire-domain  basis  functions  on  a  given  body  of  revolution. 

Third,  actual  range  measurements  could  be  made  under  the  best  of  controlled  con¬ 
ditions  following  computation  of  the  scattering  by  this  entire-domain  method  for  various 
bodies  of  revolution.  Quantification  of  the  actual  number  of  entire-domain  modes  required 
to  achieve  given  tolerances  could  be  done. 

Fourth,  an  error  analysis  study  could  be  performed  based  on  the  fact  that  the  entire- 
domain  functions  contain  frequency  content  information,  and  thus,  the  expansion  coeffi¬ 
cients  are  related  to  the  frequency  content  of  the  unknown  current.  The  study  would  show 
the  relationship  between  the  error  and  the  number  of  entire-domain  basis  functions  used 
(something  that  cannot  be  done  with  subdomain  basis  functions) . 


6-3 


Appendix  A.  The  Computer  Code 

This  appendix  contains  MATLAB  5  computer  code  for  computing  the  scattering  cross- 
section  from  a  body  of  revolution.  The  function  GEO  defines  the  geometry  for  a  sphere 
and  can  be  rewritten  to  return  similar  parameters  for  another  geometry.  The  code  is 
initiated  by  a  call  to  MAINSPHERE  with  the  appropriate  inputs  as  outlined  in  the  code. 
The  other  functions  are  called  by  MAINSPHERE. 


xxxxxxxxx 


:xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


function  [sdout]  =  mainsphere(pl,  p2,  p3,  p4,  p5,  p6) 

*/, usage:  [sdout]  =  mainsphere(a_samples,  sdmax,  fmin,  fmax,  theta,  ’ed’  or 


txxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


’sd’) 


'/.'X  This  in  the  main  program  for  the  entire  (or  sub)  domain  analysis  %’/, 


axxxxxxxxxxxxxx  mm 


VXXXXXXXXXXXXXXV,J 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XXXXX  INITIALIZE  VARIABLES  UW, 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


global  k  a  skosh  '/,  geo  and  other  constants 
global  fidd  '/,  output  control,  l=display 
global  n  ntest  nbasis  '/,  mode  indices 
global  tlo  thi  tplo  tphi  '/,  integration  limits 
global  theta_t  */,  monostatic  angle 
global  T  '/,  basis  polynomial  table 

global  bpx  bpy  wfxy  '/,  2d  gaussian  roots  and  weights  (for  v) 
global  bpx3  bpy3  bpz3  wfxyz  '/,  3d  gaussian  roots  and  weights 


(for  z) 


if  p6==,sd> 

T  =  zeros (25, 25) ; 

T(:,25)  =  [1];  '/,  pulse  basis  polynomial  coefficient 
else 

T  =  maket;  X  make  the  Chebyshev  basis  lookup  table 
end; 

fidd  =  1;  '/,  display 

lambda  =  1;  '/,  wavelength  (1  =  normalized) 
k  =  2*pi/lambda;  '/,  wave  number 


xxxxxxxxxxxxxxxxxxxxxxx 

van  CHANGEABLES  XXXXX 
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xxxxxxxxxxxxxxxxxxxxxxx 

a_samples  =  pi;  7  sphere  radii  vector 

sdmax  =  p2;  7,  initial  max  sub-domain  mode 

fmin  =  p3;  7*  min  fourier  mode  (set  ~0  for  partial  computation) 

fmax  =  p4;  7#  max  fourier  mode 

theta  =  p5;  7*  increments  in  theta 

zpoints  =  [30  30  30]  ;  7.  #  of  integration  points  for  z-matrix 

vpoints  =  [30  30]  ;  7.  #  of  integration  points  for  v-matrix 

7.fname  (below)  7#  save  filename. mat 
7*fdir  (below)  7  save  directory 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

77.7.7.7.  ENTIRE  DOMAIN,  MOMENT  METHOD  77.7.7.7. 
77777777.7.7.777.7777.77.77.77.7.7.77 


for  aindex=l :length(a_samples) 

a  «  a_samples(aindex)/k;  7  sphere  radius 
tmax  =  pi*a;  7  total  arclength  of  generating  curve 
skosh  =  .01; 

fname  =  [’sp’  num2str(f loor(k*a*10))  p6  num2str (sdmax)  ’f’  num2str(fmax)] ; 
fdir  =  ’~/outf iles/’ ; 

eval([’save  * ,  fdir,  fname,  J  a  sdmax  fmin  fmax  theta  zpoints  vpoints;’]); 

if  p6==’sd’7  make  array  of  integration  limits 
limits  =  makelimits(tmax,  sdmax); 
else 

tlo  =  0; 
thi  =  pi*a; 
tplo  =  tlo; 
tphi  =  thi; 
end; 

[bpx  bpy  wfxy]  =  grule2d(vpoints(l) ,  vpoints(2)); 

[bpx3  bpy3  bpz3  wfxyz]  =  grule3d(zpoints(l) ,  zpoints(2),  zpoints(3)); 

for  nfourier  =  fmin: fmax  7  fourier  mode  loop 

makezv;  7  make  Z,  V,  R,  I 

7  a  test  In  for  convergence  here 

7  Store  variables 

ns  =  num2str (nfourier) ; 

eval([’Rt’  ns  ’  =  transpose (Rt) ;’]) ; 
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eval([’Rp’  ns  ’  =  transpose(Rp) ;  ’]  )  ; 
eval([’I’  ns  ’=In;’]); 
eval ( [ ’ z ’  ns  * =z ; ’ ] ) ; 
eval([’v’  ns  ’=v;’]); 

eval([’save  ’ ,  fdir,  fname,  ’  Rt’  ns  9  Rp’  ns  9  I * ... 
ns  ’  z9  ns  ’  v’  ns  9  -append;’]); 

end;  X  for  nfourier 

if  fmin==0 

comprcs;  X  compute  res 

eval ([’save  9 ,  fdir,  fname,  ’  sigma_lam2  sigma_db  -append;’]); 
plotrcs;  X  plot  res 
sdout  =  sigma_db; 
else 

disp(’Partial  series  of  f-modes  computed,  no  sigma  calculated.’); 
sdout  =  9999; 
end;  X  if 

end;  X  for  aindex 


xxxxxxxxxxxxxxxxxxxxxxxx 


X  Procedure:  makezv.m 

X  Makes  the  z  and  v  matricies  and  computes  I 
X  Also  makes  the  Rt  and  Rp  matricies 


t_0  =  eputime; 


n  =  nfourier; 

A  =  zeros (sdmax+1) ; 
B  =  A; 

C  *  A; 

D  =  A; 


E  =  zeros (sdmax+1,  length (theta) ) ; 
F  =  E; 


XXXXXXXXXXX 

X  Compute  for  initial  number  of  sd  modes 
for  ntest  =  Orsdmax 

for  nbasis  =  0:sdmax 
if  p6==’sd’ 

tlo  =  limits(ntest+l,l) ; 
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thi  =  limits (ntest+1, 2) ; 
tplo  =  limits (nbasis+1, 1) ; 
tphi  =  limits(nbasis+l,2) ; 

end;  7.  if 

fprintf (f idd,  ’\n - ztp - \n’); 

tc=cputime;  B(ntest+1, nbasis+1)  =  ztp2;  t_cpu=cputime-tc; 
fprintf (f idd, ’\nElapsed  cpu  time:  7,f  seconds. \n’  ,t_cpu) ; 

if  nbasis  <=  ntest 

fprintf  (f  idd,  ’  \n - ztt - \n’)  ; 

tc=cputime;  A(ntest+1, nbasis+1)  =  ztt2;  t_cpu=cputime-tc; 
fprintf  (f  idd, ’\nElapsed  cpu  time:  7,f  seconds. \n ’ ,t_cpu) ; 

fprintf  (f  idd,  ’\n - zpp - \n’)  ; 

tc=cputime;  D(ntest+1, nbasis+1)  =  zpp2;  t_cpu=cputime-tc; 
fprintf  (f  idd,  J\nElapsed  cpu  time:  7,f  seconds. \n’ ,t_cpu) ; 

end;  7,  if 

end;  7,  for  nbasis 

for  th  =  l:length(theta) ; 

theta_t  =  theta (th) ; 

f  pr  int f  (  f  idd ,  ’  \n - vt - \n  ’ )  ; 

tc=cputime;  E(ntest+1,  th)  =  vt;  t_cpu=cputime-tc; 
fprintf  (f  idd, ’\nElapsed  cpu  time:  7.f  seconds  An’  ,t_cpu)  ; 

fprintf  (f  idd,  ’\n - vp - \n’)  ; 

tc=cputime;  F(ntest+1,  th)  =  vp;  t_cpu=cputime-tc; 
fprintf  (f  idd,  *\nElapsed  cpu  time:  7.f  seconds. \n’ ,t_cpu)  ; 

n  =  -n; 

Rt (ntest+1,  th)  =  vt; 

Rp (ntest+1,  th)  =  vp; 

n  =  -n;  7.  restore 

end;  7.  for  th 
end;  7.  for  ntest 

A  =  tril(A)  +  transpose (tril (A,  -1));  7.  exploit  symmetries 
D  =  tril(D)  +  transpose (tril (D,  -1)); 

C  =  -transpose (B) ; 

z  =  [A ,  B ;  C ,  D]  ; 
v  =  [E;  F]  ; 
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fprintf  (fidd,  *\n - Invert - \n’) ; 

tc= cputime;  In  =  inv(z)  *  v;  t_cpu=cputime-tc; 

fprintf (fidd, ’XnElapsed  cpu  time:  #/,f  seconds  An’  ,t_cpu)  ; 

t_tot=cputime  -  t_0; 

disp([* Total  runtime:  *  num2str(t_tot)] ) ; 
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X  Procedure:  comprcs.m 


fprintf  (fidd,  *\n - Compute  RCS - \n*)  ; 

for  th  =  1: length (theta) 
theta_t  =  theta(th); 
nsum  =  0; 
for  n  =  l:fmax 

fprintf  (fidd,  *\nAdd  n=,/,d\n>,  n)  ; 

ns  =  num2str(n) ; 

eval([’In  =  I9  ns  *  (:,  th);*]); 

len  =  length  (In); 

mid  =  len/2; 

It  =  In(l:mid) ; 

Ip  =  In(mid+l:len) ; 

nsum  =  nsum  +  eval([*Rt*  ns  *  (th,  :)  *  It  +  Rp*  ns  *  (th,  :)  *  Ip*]) 
end;  '/,  for  n 


In  =  10 ( : ,th) ; 
len  =  length(In); 
mid  =  len/2; 

ItO  =  In(l  :  mid); 

IpO  =  In(mid+1  :  len) ; 

sigma_lam2(th)  =  l/(4*pi~3)  *  (  abs(.5  *  Rt0(th,  :)  *  ItO  +  ... 
.5  *  Rp0(th,  :)  *  IpO  +  nsum)  )~2; 

end;  X  for  th 

if  exist( ,aindex*)==l 

sigma.db(aindex)  =  10*logl0(sigma_lam2) ; 
else 

sigma_db=10*logl0(sigma_lam2)  ; 
end; 
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*/,  Procedure:  plotrcs*m 


figure; 

plot(theta*180/pi,  sigma.db); 

titled ’Scattering  Cross-Section’ , [’Date:  ’ ,  date,  ’  Datafile 
fdir,  ’  ’,  fname]>); 

sdata  =  num2str([sigma_db(l:10);  sigma_db(ll :20)] ) ; 
xlabel ({’Theta  (deg)’,  sdata}); 
ylabel(’Sigma/lambda~2  (dBsm)’); 
grid  on; 
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function  [A]  =  makelimits(tmax,  sdmax) 

X  makes  an  n  x  2  array  of  integration  limits 
X  for  the  subdomain  analysis 


X  as  in  entire  domain,  the  modes  include  the  O-mode 

delta  =  tmax  /  (sdmax  +  1); 

lo(l)  =  0; 

hi(l)  =  delta; 

for  index  =  2:sdmax+l 

lo( index)  =  hi ( index- 1); 
hi(index)  =  lo(index)  +  delta; 
end; 


A  =  [lo\  hi’]; 


XXXXXXXXXXXXXXXXXXXXXXXX 


function  [A]  =  maket 

X  returns  the  vector  of  coefficients 

X  of  the  Chebyshev  polynomial  of  mode  "mode11 


hi=25; 

t=zeros(hi) ; 


t  (1  ,hi)  =  [1]; 
t(2,hi-i:hi)  =  [10]; 
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for  nn=3 :25 

s  =  conv([2  0] ,t(nn-l, :)) ; 
t(nn,l:hi)  =  s(2:hi+l)  -  t(nn-2,:); 
end;  7,  for 

A  =  t; 
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function  [bpx,bpy,wfxy]  =  grule2d  (nquadx,nquady) 


[bpxv,wfxv] =grule(nquadx) ; 
[bpyv,wfyv]=grule(nquady) ; 
[bpx,bpy]=meshgrid(bpxv,bpyv) ; 
[wfx.wfy] =meshgrid(wfxv,wfyv) ; 
wf xy=wf x . *wf y ; 


xxxxxxxxxxxxxxxxxxxxxxxx 


function  [bpx,bpy,bpz,wfxyz] 


grule3d  (nquadx , nquady , nquadz ) 


[bpxv , wf xv] =grule (nquadx) ; 

[bpyv , wf yv] =grule (nquady) ; 

[bpzv.wfzv] =grule (nquadz) ; 

[bpx ,bpy ,bpz] =meshgr id (bpxv , bpyv ,bpzv) ; 
[wfx,wfy,wfz]=meshgrid(wfxv,wfyv,wfzv) ; 
wf xyz=wf x . *wf y . *wf z ; 


xxxxxxxxxxxxxxxxxxxxxxxx 


function  [bp,wf]=grule(n) 

X  This  function  computes  Gauss  base  points  and  weight  factors 
X  using  the  algorithm  given  by  Davis  and  Rabinowitz  in  ’Methods 
X.  of  Numerical  Integration’,  page  365,  Academic  Press,  1975. 

X  by  Howard  Wilson  of  the  University  of  Alabama,  1990 
X  included  in  the  Numerical  Integration  Toolbox* 

X  *duplicated  here  for  convienence  of  the  reader 
bp=zeros(n,l) ;  wf=bp;  iter=2;  m=f ix((n+l)/2) ;  el=n*(n+l); 
mm=4*m-l;  t=(pi/(4*n+2))*(3:4:mm) ;  nn=(l-(l-l/n)/(8*n*n)) ; 
xo=nn*cos(t) ; 
for  j=l liter 

pkml=l ;  pk=xo; 
for  k=2:n 
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t 1 =xo . *pk ;  pkp 1 =t 1 -pkm 1 - ( t 1 -pkml ) /k+t 1 ; 
pkml=pk;  pk=pkpl; 
end 

den=l.-xo.*xo;  dl=n*(pkml-xo.*pk) ;  dpn=dl./den; 
d2pn=(2.*xo.*dpn-el.*pk) ./den; 
d3pn=(4*xo.*d2pn+(2-el) .*dpn) ./den; 
d4pn=(6*xo.*d3pn+(6-el) .*d2pn) ./den; 
u=pk . / dpn ;  v=d2pn . / dpn ; 

h=-u.*(l+( ,5*u) .*(v+u.*(v.*v-u.*d3pn./(3*dpn)))) ; 
p=pk+h . * (dpn+ ( . 5*h) . * (d2pn+ (h/3) . * (d3pn+ . 25*h . *d4pn) ) ) ; 
dp=dpn+h . * (d2pn+ ( . 5*h) . * (d3pn+h . *d4pn/3) ) ; 
h=h-p./dp;  xo=xo+h; 

end 

bp=-xo-h ; 

fx=dl-h.*el ,*(pk+(h/2) ,*(dpn+(h/3) .*(... 

d2pn+(h/4) .*(d3pn+( .2*h) .*d4pn)))) ; 
wf =2* (1-bp. “2) ./(fx.*fx) ; 
if  (m+m)  >  n,  bp(m)=0;  end 
if  ~((m+m)  ==  n) ,  m=m-l ;  end 

j j=l :m;  nlj=(n+l-j j) ;  bp(nlj )=-bp(j j ) ;  wf (nlj)=wf (j j) ; 

7,  end 


mxxmxxxxmmmm 


function  vol  =  gquad2d(fun,xlow,xhigh,ylow ,yhigh,bpx,bpy ,wfxy) 

7, based  on  a  routine  by  Bryce  Gardner,  Purdue  University,  Spring  1993 


7.Map  to  x 
qx= (xhigh-xlow) /2 ; 
px= (xhigh+xlow) /2 ; 
x=qx*bpx+px ; 

7.Map  to  y 
qy= (yhigh-ylow) /2 ; 
py= (yhigh+ylow) /2 ; 
y=qy*bpy+py; 

fv  =  f eval(fun,x,y) ; 

vol  =  sum(sum(wfxy.*fv))*qx*qy; 

f  printf  ( 1 ,  *  The  7.d  point  integration  of  * ’Xs”  gives\n  7.20. 15f  ... 
+  i7.20.15f\n’ ,  lengthCbpx),  fun,  real(vol),  imag(vol)); 
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function  vol  =  gquad3d(fun,xlow,xhigh,ylow,yhigh,zlow,zhigh,bpx,bpy ,bpz,wfxyz) 
'/based  on  a  routine  by  Bryce  Gardner,  Purdue  University,  Spring  1993 


'/Map  to  x 
qx= (xhigh-xlow) / 2 ; 
px= (xhigh+xlow) / 2 ; 
x=qx*bpx+px ; 


'/Map  to  y 
qy= (yhigh-ylow) /2 ; 
py= (yhigh+ylow) /2 ; 
y=qy*bpy+py; 

'/Map  to  z 

qz=  (zhigh-zlow)  /2 ; 
pz= (zhigh+zlow) /2 ; 
z=qz*bpz+pz ; 

fv  =  f eval(fun,x,y,z) ; 

vol  =  sum(sum(sum(wfxyz.*fv)))*qx*qy*qz; 

fprintf  (1, 'The  '/,d  point  integration  of  9  9%s9*  gives\n'/,20*  15f 
+  i'/20.15f\n'  ,length(bpx)  , fun, real  (vol)  ,imag(vol))  ; 


xxxxxx  mmm 
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function  [A]  =  basis (mode,  tp) 

'/  Returns  the  value  of  the  basis  function 
X  of  mode  'mode'  evaluated  at  tp. 
global  T  tphi 


s  =  2  *  tp  ./  tphi  -  1;  */  scales  domain  to  +-1 


'/  workaround  for  3d  (nul,  nu2  are  not  needed) 
Al=polyval(T(mode+l, :),s(:,l,l)); 

[nul,  A,  nu2]=meshgrid(Al,  Al,  Al) ; 


function  [A]  =  derivtp(mode,  tp) 


'/  Computes  the  partial  derivative  with  respect  to  tp 
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'/,  of  rho(tp)  *  basis  (tp) 

•/.  assumes  sphere 
global  a 

A  =  cos(tp/a)  .*  basis (mode,  tp)  +  a*sin(tp/a)  .*  dbasis(mode,  tp) ; 


mmv/x/T 


function  [A]  =  dbasis(mode,  tp) 

X  returns  the  value  of  the  derivative 
*/,  of  the  basis  function  @  mode,  tp 
global  T  tphi 


s  =  2  *  tp  ./  tphi  -  1;  '/,  scales  domain  to  +-1 

'/,  workaround  for  3d  (nul,  nu2  are  not  needed) 

A1  =  polyval(polyder(T(mode+l , : )) ,s( : , 1 ,1) ) ; 
[nul,  A,  nu2]=meshgrid(Al,  Al,  Al); 


xxxxxxxxxxxxxxxxxxxxxxxxx 


function  [rho,  zee,  vee]  =  geo(t) 

*/,  Computes  the  geometric  parameters  for  a 
X  ***  SPHERE  *** 

X  all  positive  z 
global  a  skosh 


next_t  =  t  +  ,00001; 

rho  =  a  *  sin(t/a); 

rho^next  =  a  *  sin(next_t/a) ; 

zee  =  -a  *  cos(t/a)  +  a; 

zee_next  =  -a  *  cos (next _t /a)  +  a; 

drho  =  rho_next  -  rho; 

dzee  =  zee_next  -  zee; 

hyp  =  sqrt(drho.~2  +  dzee, "2); 

vee  =  as in (drho. /hyp ) ; 


xxxxxxxxxxxxxxxxxxxxxxxxxxxx 


function  [A]  =  test (mode,  t) 

X  Returns  the  value  of  the  testing  function 
X  of  mode  ’mode’  evaluated  at  t. 
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'/,  Galerkin  (basis  =  test) 
global  T  thi 

s  =  2  *  t  ./  thi  -  1;  7.  scales  domain  to  +-1 

if  length(size(s))  ==  3 

7.  workaround  for  3d  (nul,  nu2  are  not  needed) 
Ai=polyval(T(mode+l, :),s(l, : ,1)); 

[A,  nul,  nu2]=meshgrid(Al,  Al,  Al) ; 
else 

A=polyval(T(mode+l , :) ,s) ; 
end; 


f miction  [A]  =  derivt (mode ,  t) 

'/,  computes  the  partial  derivative  with  respect 
'/,  to  t  of  rho(t)  *  test(t). 

'/,  assumes  sphere 
global  a 


A  =  cos(t/a)  .  *  test (mode,  t)  +  a*sin(t/a) .*dtest(mode,  t) ; 


mm 


function  [A]  =  dtest(mode,  t) 

'/,  returns  the  value  of  the  derivative 
7,  of  the  test  function  @  mode,  t 
'/,  galerkin:  test  =  basis 
global  T  thi 


s  =  2  *  t  ./  thi  -  1;  */,  scales  domain  to  +-1 


'/,  workaround  for  3d  (nul,  nu2  are  not  needed) 
Al  =  polyval(polyder(T(mode+l, :)) ,s(l, : ,1)) ; 
[A,  nul,  nu2]=meshgrid(Al,  Al,  Al); 


mmmmmmmmram 


function  [A]  =  ztt2 

'/,  Computes  an  impedance  matrix  element  of  ztt 
global  tplo  tphi  tlo  thi 
global  bpx3  bpy3  bpz3  wfxyz 
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A  =  gquad3d(’zttarg2’ ,  tlo,  thi,  tplo,  tphi,  0,  pi,  bpx3,  bpy3,  bpz3,  wfxyz); 


function  [A]  =  zttarg2(t,tp,phip) 
'/,  Argument  for  ztt  integration 
global  k  n  skosh  ntest  nbasis 


[rho_t,  zee_t,  vee_t]  =  geo(t); 
[rho_tp,  zee_tp,  vee_tp]  =  geo(tp) ; 


R  =  sqrt(  skosh. “2  +  (rho_t-rho_tp) . ~2  +  (zee_t-zee_tp) .“2  +  ... 

4*rho_t.*rho_tp.*(sin(phip/2)) ."2  ) ; 

arg4  *  exp(-j*k*R) ./(k*R)  .*  cos(n*phip); 

arg5  =  exp(-j*k*R) ./(k*R)  .*  cos(phip)  .*  cos(n*phip); 

A1  =  j*k“2*rho_t.*test(ntest,t) .*rho_tp.*basis (nbasis, tp) .*  ... 

(  arg5.*sin(vee_t) . *sin(vee_tp)  +  arg4. *cos(vee_t) . *cos(vee_tp)  ); 
A2  =  -j .*derivt (ntest, t) .*derivtp(nbasis,tp) .*arg4; 

A  =  A1  +  A2; 


mmmnm 


mm 


function  [A]  =  ztp2 

'/,  Computes  an  impedance  matrix  element  of  ztp 
global  tplo  tphi  tlo  thi 
global  bpx3  bpy3  bpz3  wfxyz 


A  =  gquad3d(’ztparg2’ ,  tlo,  thi,  tplo,  tphi,  0,  pi,  bpx3,  bpy3,  bpz3,  wfxyz); 


•/.mm 


function  [A]  =  ztparg2(t,tp,phip) 
'/,  Argument  for  ztp  integration 
global  k  n  skosh  ntest  nbasis 


[rho_t,  zee_t,  vee_t]  =  geo(t); 
[rho_tp,  zee_tp,  vee_tp]  =  geo(tp) ; 


R  =  sqrtC  skosh. "2  +  (rho_t-rho_tp) .“2  +  (zee_t-zee_tp) ."2  +  ... 
4*rho_t.*rho_tp.*(sin(phip/2)) . ~2  ) ; 
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arg4  =  exp(-j*k*R) ./(k*R)  .  *  cos(n*phip); 

arg6  =  exp(-j*k*R) ./(k*R)  .*  sin(phip)  .*  sin(n*phip) ; 

A1  =  rho_tp.*basis (nbasis, tp) . *k~2.*rho_t . *test (ntest, t) .*arg6.*sin(vee_t) ; 
A2  ®  rho_tp.*basis(nbasis,tp) .*  n./rho_tp  .*derivt (ntest, t) .*arg4; 

A  =  A1  +  A2; 


xxxxxxxxxxxxxxxxxxx 


'XXXXXX 


function  [A]  =  zpt2 

X  Computes  an  impedance  matrix  element  of  zpt 
global  tplo  tphi  tlo  thi 
global  bpx3  bpy3  bpz3  wfxyz 


A  =  gquad3d(,zptarg2> ,  tlo,  thi,  tplo,  tphi,  0,  pi,  bpx3,  bpy3,  bpz3,  wfxyz); 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


function  [A]  =  zptarg2(t ,tp,phip) 
*/,  Argument  for  zpt  integration 
global  k  n  skosh  ntest  nbasis 


[rho_t,  zee_t,  vee_t]  =  geo(t); 

[rho.tp,  zee_tp,  vee_tp]  =  geo(tp); 

R  =  sqrt(  skosh. "2  +  (rho_t-rho_tp) . ~2  +  (zee_t-zee_tp) . ~2  +  ... 

4*rho_t . *rho_tp.*(sin(phip/2)) .  ~2  ); 

arg4  =  exp (- j  *k*R) . / (k*R)  .*  cos(n*phip); 

arg6  =  exp(-j*k*R) ./(k*R)  . *  sin(phip)  .*  sin(n*phip); 

A1  =  -rho_t.*test (ntest,  t) . *k~2.*rho_tp. *basis (nbasis, tp) .*arg6 . *sin(vee_tp) ; 
A2  =  ~rho_t.*test (ntest,  t).*  n./rho_t  . *derivtp (nbasis, tp) . *arg4; 

A  =  A1  +  A2; 


function  [A]  =  zpp2 

*/,  Computes  an  impedance  matrix  element  of  zpp 
global  tplo  tphi  tlo  thi 
global  bpx3  bpy3  bpz3  wfxyz 


A  =  gquad3d( ,zpparg2> ,  tlo,  thi,  tplo,  tphi,  0,  pi,  bpx3,  bpy3,  bpz3,  wfxyz); 
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function  [A]  =  zpparg2(t ,tp,phip) 
l  Argument  for  zpp  integration 
global  k  n  skosh  ntest  nbasis 


[rho_t,  zee_t,  vee_t]  =  geo(t); 
[rho_tp,  zee_tp,  vee_tp]  =  geo(tp); 


R  =  sqrt(  skosh. "2  +  (rho_t-rho_tp) ."2  +  (zee_t-zee_tp) .~2  +  ... 

4*rho__t .  *rho_tp .  * ( sin (phip/2) ) .  "2  ) ; 

arg4  =  exp(-j*k*R) ./(k*R)  .*  cos(n*phip); 

arg5  =  exp(-j*k*R) ./(k*R)  .  *  cos(phip)  .*  cos(n*phip); 

A  =  j *rho_t .*test (ntest, t) .*rho_tp.*basis(nbasis,tp) .  *  ... 

(  k~2.*arg5  -  n~2./(rho_t . *rho_tp)  .*  arg4  ); 


function  [A]  =  vt 

l  Computes  a  voltage  matrix  element  of  vt 
global  tlo  thi 
global  bpx  bpy  wfxy 


A  »  gquad2d(,vtarg> ,  tlo,  thi,  0,  2*pi,  bpx,  bpy,  wfxy); 


mr 


function  [A]  =  vtarg(t,phi) 

7,  argument  for  vt  integral, 
global  k  n  ntest  theta_t 


[rho,  zee,  vee]  =  geo(t); 
dphi  =  sin(vee)  .*  sin(phi) ; 

V.dtheta  =  sin(vee(t))  .  *  cos  (phi); 

minus  _kr  =  k  *  zee  .  *  cos(theta_t)  +  k  *  rho  .  *  sin(theta_t)  .*  cos  (phi); 
preA  =  k  *  rho  .*  test (ntest ,t)  .*  exp(j .*(minus_kr  -  n.*phi)); 

A  =  preA  .*  dphi; 
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function  [A]  =  vp 

'/,  Computes  a  voltage  matrix  element  of  vp 
global  tlo  thi 
global  bpx  bpy  wfxy 


A  =  gquad2d( ’ vparg ’ ,  tlo,  thi,  0,  2*pi,  bpx,  bpy,  wfxy); 


mtmtnm 


mi 


function  [A]  =  vparg (t, phi) 

'/,  argument  for  vp  integral, 
global  k  n  ntest  theta_t 


[rho,  zee,  vee]  =  geo(t); 
dphi  =  cos (phi); 

‘/.dtheta  =  -sin(phi); 

minus.kr  =  k  *  zee  .  *  cos(theta_t)  +  k  *  rho  .*  sin(theta.t)  .*  cos (phi) 
preA  =  k  *  rho  .*  test(ntest.t)  .*  exp(j*(minus_kr  -  n.*phi)); 

A  =  preA  .*  dphi; 
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Appendix  B.  Current  Densities  by  Mode 

This  appendix  contains  the  current  density  components  for  all  modes  of  the  ka  =  4.0 
sphere  in  both  the  t  and  0  directions  for  E™c  incident  at  0  =  90°.  The  total  current 
density  over  accumulated  over  the  Tq-T$  Chebyshev  modes  and  the  n  =  0-n  =  5  is  also 
included. 


^Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=0,  E^c 
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^-Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=3,  E^nc 
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(^-Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=5,  E^nc 
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t-Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=2,  E^nc 
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t-Component  of  Current  Density  for  0=90  deg  and  ka  =  4,  n=5,  E™ 
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Total  current  density  over  all  ed/f  modes,  0=90  deg,  ka=4 
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